Cluster-CPU-GPU heterogeneous parallel method for solving response of large-scale anisotropic structure

Abstract

The invention discloses a cluster-CPU-GPU heterogeneous parallel method for solving response of a large-scale anisotropic material composition structure, and belongs to the technical field of high-performance calculation simulation. The method comprises the following steps: firstly, establishing a macroscopic structure and a microstructure, dividing a finite element mesh for a macroscopic model, dividing a calculation region by using a region decomposition algorithm, and realizing macroscopic model unit stiffness matrix and displacement calculation through OpenMP acceleration; secondly, reconstructing the model at the Gaussian point of the macroscopic model, partitioning by means of the MPI technology, calculating the displacement of the microscopic model by using multiple threads in each process, and accelerating the calculation of a stiffness matrix of a microscopic model unit by using OpenMP to realize cluster MPI parallel calculation; during post-processing, visualization is realized by utilizing GPU multi-core advantage calculation macroscopic model stress and strain, and cross-scale heterogeneous parallel simulation is achieved. The method is combined with related algorithms, can efficiently utilize resources to solve the analog simulation problem, and is wide in application prospect.

Description

Technical Field

[0001] The present invention belongs to the technical field of high-performance finite element simulation, and specifically is a cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural responses. Technical Background

[0002] Anisotropic materials have the characteristic that their physical properties show certain differences in different directions. Anisotropic materials are mainly divided into crystals, composite materials, fiber materials, etc. The research on these materials needs to consider their unique properties. Anisotropic materials are currently also widely used in medical, aerospace and other fields.

[0003] For a long time, the research on the macroscopic properties of anisotropic materials has been mainly based on experimental trial and error. However, this method is not only costly, has a large chance and a long test cycle, but also a huge waste of manpower and material resources. At present, with the rapid development of computer technology, the use of high-performance computers to predict anisotropic materials has become a hot topic of research. Compared with traditional experimental methods, the use of computers for numerical simulation calculations can not only design the microstructure of materials more conveniently, but also be simpler and more repeatable than experimental methods, and can more quickly derive the mechanical properties of materials, so as to facilitate the next step of optimization and improvement. In order to obtain the mechanical properties of materials more accurately, a larger-scale model is required, but large-scale solutions often affect the accuracy and correctness of the solutions. At present, the main frequency of computers has increased less, and the number of cores has increased more. Only by making full use of computing resources can higher-precision and efficient solutions be achieved.

[0004] Therefore, a cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural responses based on regional decomposition is established, and the mature finite element numerical calculation theory is applied. The advantages of different parallel methods are fully utilized to accelerate the calculation. It can more effectively use computing resources to achieve rapid solutions to large-scale models. At the same time, the cross-scale method 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 aspects, and promote the optimization of mechanical properties of materials and structural design. Summary of the invention

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

[0006] The present invention provides a solution to the above technical problem: a cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response, comprising the following steps:

[0007] S1. Establish any large-scale macro-geometric model and micro-geometric model, discretize the macro-geometric model and micro-geometric model to be analyzed, divide the grid, and generate the required grid data;

[0008] S2. Set boundary conditions, constraints, and material parameters for the macro model, and use CPU multi-threading technology (OpenMP) to parallelly calculate the unit stiffness matrix of the macro model.

[0009] S3. Use the regional decomposition algorithm to partition the model according to the set value, divide the grid model to be solved into multiple calculation sub-regions and make each partition carry all its own parameters. Iterate in each partition according to the partition information and solve the displacement field of the complete model. During the solution process, use the open source library STXXL to place the data in the memory on the hard disk, so as to reduce memory usage and quickly read the hard disk data.

[0010] S4. According to the calculated displacement results, the strain and stress at each Gauss point of each unit of the macro model are calculated;

[0011] S5. Use the multi-process parallel technology (MPI) of the cluster to evenly distribute the units of the macro model 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 Gauss point of each unit of the macro model under each process, use the micro model as the microstructure of the macro model, and use the strain at the Gauss point to calculate the displacement boundary condition of the microstructure at the point;

[0013] S7, using CPU multi-threading technology (OpenMP) in each process to parallelly calculate the unit stiffness matrix of each unit on the microstructure, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each unit of the micromodel according to 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) according to the stress of each Gaussian point on the micro-model as the average stress of a single Gaussian point of a single unit of the macro-model, and then the internal force of the unit of the macro-model is calculated according to the average stress of four Gaussian points of a single unit of the macro-model;

[0015] S9, each process completes the internal force calculation of all units in the process, process 0 collects the calculation results of the internal forces of all units of the macro model to obtain the internal force vector of the macro model, and calculates the displacement of the macro model according to the internal force vector;

[0016] S10, iterating the displacement judgment residual obtained by using the displacement calculated by the displacement boundary of the macro model and the internal force vector calculation, and obtaining the displacement result of the macro model after judging convergence;

[0017] S11. Calculate the stress and strain of the macro model based on the displacement results through GPU, quickly implement post-processing of the calculation, and realize the cross-scale heterogeneous parallel algorithm of the anisotropic material composition structure.

[0018] The specific implementation method of step S1 is: establish a three-dimensional tetrahedral basic grid unit with a reasonable density according to the complexity of the macroscopic scale and the microscopic scale. The model saves the number of model units, the nodes contained in each unit, and the number of nodes and the three-dimensional coordinates of all nodes.

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

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

[0021] Set the material parameters of the macro model and micro model respectively. The macro model is regarded as an isotropic material, and the Poisson's ratio and elastic modulus are set according to the requirements. The micro model is an anisotropic material, and the elastic modulus and Poisson's ratio in each direction need to be set. For the macro model, it is also necessary to set boundary conditions and constraints to determine the location of the displacement boundary and fixed constraints.

[0022] S22, calculation of element stiffness matrix;

[0023] According to the set macro model material parameters, the macro model unit stiffness matrix is ​​calculated. First, the shape function of the three-dimensional tetrahedral unit is calculated. For the tetrahedral unit used, its volume coordinate is defined as;

[0024]

[0025]

[0026] The shape function in the parametric coordinate system is equal to the corresponding volume coordinate. The four shape functions corresponding to the four nodes are:

[0027]

[0028] The partial derivative of the shape function in each direction of the local coordinate system is then multiplied by the node coordinates to obtain the Jacobian matrix. The method for solving the Jacobian matrix is ​​as follows;

[0029]

[0030] Multiply the inverse matrix of the Jacobi matrix by the partial derivative of the shape function to get the B matrix;

[0031]

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

[0033]

[0034] Among them, E and v are the elastic modulus and Poisson’s ratio of the macro model, respectively.

[0035] The expression of the unit stiffness matrix is ​​obtained by deducing the energy principle. Since the expression of the stiffness matrix is ​​an integral operation, the calculation of the unit stiffness matrix is ​​realized by introducing numerical integration. The results of hammer numerical integration are as follows;

[0036]

[0037] S23, OpenMP parallel computing;

[0038] Since the calculation of the element stiffness matrix of each element is independent, and the calculation of the element stiffness matrix is ​​only performed between a single loop, and involves fewer variables, the OpenMP multi-threading technology is used to perform parallel calculation of the element stiffness matrix. The elements of the macro model are divided into blocks, and each thread calculates the element stiffness matrix of the element assigned to each thread.

[0039] At this point, the setting of boundary conditions, constraints, and material parameters for the macro model has been completed, and the unit stiffness matrix of the macro model is calculated in parallel using CPU multi-threading technology (0penMP).

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

[0041] S31, model domain decomposition;

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

[0043] After the partitioning is completed, 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 the renumbering will be saved. At this time, the original finite element model is divided into N blocks of regions Ω0~Ω N-1

[0044] S32, partition stiffness matrix calculation;

[0045] After completing the regional division of the model, the partition stiffness matrix calculation of each partition is integrated according to the partition situation. According to the partition unit information in S31, the unit stiffness matrix of each partition containing unit is extracted, and the partition stiffness matrix is ​​constructed according to the printed node degree of freedom encoding of the partition re-edited node information.

[0046] S33, iteratively solving the displacement field;

[0047] First, the global Schwarz solution iteration in the regional decomposition algorithm is:

[0048]

[0049] where ω is an arbitrary function that maps the domain Ω to R, and u n is an approximate solution to the Poisson problem. This formula is used to solve the subproblems in all subdomains of the regional algorithm at the same time, and its convergence is proved using the maximum principle. On this basis, the residual is introduced as follows;

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

[0051] The restricted additive Schwarz algorithm (RAS) on the continuous level can be defined.

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

[0053] KU=F

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

[0055]

[0056] Since the RAS algorithm is defined based on the continuous level, and 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, let D be the diagonal matrix of K, then;

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

[0058] Define the residual at the discrete level as r n =F-KU n , the Jacobi algorithm is divided into blocks, and the block Jacobi algorithm can be defined as;

[0059]

[0060] Among them, K ii is the overall stiffness matrix of each partition, is the displacement vector of the nth iteration step of the i-th partition. At this time, R is introduced i As N to N i The restriction operator, As N to N i The transpose operator of At this time, the block Jacobi algorithm can be written as follows:

[0061]

[0062] Based on the above formula, the RAS algorithm applicable to any number of partitions at the discrete level is:

[0063]

[0064] For the block Jacobi algorithm, the matrix expansion with the partition number of 2 is:

[0065]

[0066] in, Represents the displacement of the first partition at the nth iteration step The non-overlapping part of Represents the displacement of the first partition at the nth iteration step overlapping part. Represents the residual of the first partition in the nth iteration The non-overlapping part, Represents the overlapped residuals of the first and second partitions at the nth iteration.

[0067] The final linear algebraic form obtained by simplifying the above formula is:

[0068]

[0069] in, is the local correction of the discrete level solution in the nth step,

[0070] Each partition solves its own linear equations, and after obtaining the results, the overlapping parts are multiplied by one-half and added. This form is applicable to all partition numbers. When the number of partitions is greater than 2, the regional division situation can be rewritten with reference to the above formula. It is only necessary to reflect the partition situation of the overlapping nodes in the constant term matrix. At this point, the displacement field solution of the RAS algorithm based on regional decomposition is completed.

[0071] S34, memory data storage;

[0072] For the displacement field solution process mentioned above, for models with tens of millions of units, a large amount of information needs to be stored in the memory in actual calculations, such as unit number, node coordinates and other model information, as well as the B matrix, material matrix and other finite element model information of each unit. This information will occupy a large amount of memory in the calculation. The open source library STXXL can be used to store memory data in the hard disk, thereby reducing memory usage and quickly reading hard disk data.

[0073] This open source library uses the most popular parallel disk model (PDM). In this model, file I / O is explicitly handled by the program. The model divides the M-byte data in memory into D blocks and transfers them between D hard disks to amortize the delay. Usually, the interface of the external storage algorithm assumes that it reads input from the external storage, that is, the hard disk, and writes the output to the external storage. STXXL uses the idea of ​​pipeline to equip the external memory algorithm with a new interface, allowing them to provide the output as a data stream directly to the algorithm that uses the output instead of writing it to the memory first. In addition, STXXL uses synchronous I / O calls to run asynchronously in a separate thread, that is, for each file, there is a read request queue and a write request queue and a thread. The main thread sends the request to the file queue, and the thread associated with the file executes the request in FIFO order. Using this method, the data in the memory can be stored on the hard disk, thereby reducing the memory usage.

[0074] At this point, the model is partitioned according to the set value using the regional decomposition algorithm, and the grid model is divided into multiple calculation sub-regions while each partition carries all its own parameters. According to the partition information, iterate in each partition and solve the displacement field of the complete model. During the solution process, the open source library STXXL is used to place the data in the memory on the hard disk to reduce memory usage and quickly read the hard disk data.

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

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

[0077] The stress and strain at the Gaussian point are calculated based on the three-dimensional tetrahedral unit used. The three-dimensional tetrahedral unit has four Gaussian points. 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] Among them, ε is the strain matrix at a single Gaussian point, recording the strains and shear strains in three directions, B is the B matrix at the Gaussian point, and U is the displacement matrix of each node on the unit sorted according to the unit connection order. After obtaining the strain matrix, the stress matrix is ​​calculated;

[0080] σ=Dε

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

[0082] Since then, the calculation of strain and stress at each Gaussian point of each unit in the macro model has been completed based on the calculated displacement results;

[0083] The specific implementation method of step S5 is as follows: since the displacement of the microscopic model at each Gaussian point needs to be calculated, and more calculations are involved rather than a single cycle, the MPI multi-process technology is used for parallel processing. First, the macroscopic unit is processed in blocks, and the stress, strain and displacement of each Gaussian point on each unit are distributed to each process using the MPI communication protocol.

[0084] The specific implementation method of step S6 is: reconstruct the model at each Gauss point on each unit of the macro model, and use the strain at the Gauss point to calculate the displacement boundary of the micro model. Since the finite element method first calculates the stress and strain at the Gauss point, and then uses the interpolation method to calculate the stress and strain at other positions of the unit, the micro model is regarded as a point on the macro model. The new displacement and stress and strain of the micro model at this time are calculated and then the result is regarded as the new stress at the Gauss point to achieve the purpose of cross-scale calculation.

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

[0086] S71, microscopic model boundary calculation;

[0087] After distributing the macro model to each process according to the previous step, each process needs to perform a complete micro model calculation at each Gauss point in the unit. At this time, each process contains all the model information of the micro model and the strain at each Gauss point. First, extract the boundary nodes and coordinates of the micro model, and use the strain at the Gauss point multiplied by the coordinate to obtain the displacement boundary conditions of the micro model;

[0088] Δd=εd

[0089] Where Δd is the displacement of each boundary node on the microscopic model, and d is the coordinate of each boundary node on the microscopic model. At this time, the boundary conditions of the microscopic model are obtained.

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

[0091] The same finite element theory in S22 is used to calculate the shape function and other parameters in the micro-model unit stiffness matrix. Since the macro-model is regarded as an isotropic material and the micro-model is regarded as an anisotropic material, the calculation method of the micro-model is different from that of the macro-model when calculating the material matrix D. At this time, the calculation method of the material matrix D of the micro-model is:

[0092]

[0093] According to the different forms of anisotropic materials, 36 different parameters in the matrix are determined. And the stiffness matrix of each element is calculated in the same way as in S2.

[0094] S73, OpenMP parallel computing;

[0095] According to the same method in S23, OpenMP is used to accelerate the calculation of the element stiffness matrix and realize the simultaneous scheduling of processes and threads.

[0096] S74, overall stiffness matrix assembly;

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

[0098] S75, solution of displacement, stress and strain;

[0099] According to the same method in S41, the displacement field in the microscopic model and the stress-strain value at each Gaussian point on the microscopic model are solved at this time.

[0100] At this point, the displacement, stress and strain of the micro model at a single Gauss point have been solved, and the CPU multi-threading technology (0penMP) is used in each process to parallelly calculate the unit stiffness matrix of each unit in the micro structure, assemble the overall stiffness matrix and solve the displacement field. According to the displacement results, the strain and stress at each Gauss point of each unit in the micro model are calculated.

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

[0102] S81, calculation of average stress of micro model;

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

[0104]

[0105] in, is the average stress of the model, σ ij is the stress at each Gauss point.

[0106] S82, OpenMP parallelism;

[0107] Since the average stress calculation of each unit in the microscopic model is independent of each other, the S23 method is used to accelerate the calculation of the average stress, and OpenMP is used to enable multi-threading to accelerate the loop.

[0108] S83, calculation of internal forces of micro-model units;

[0109] After obtaining the average stress of a single microscopic model, the macroscopic model is returned, that is, the stress at a single Gaussian point of the macroscopic model is equivalent to the average stress of the microscopic model at the Gaussian point, and the internal force of the macroscopic model at this time is calculated based on the average stress;

[0110]

[0111] Among them, B is the B matrix of the macro unit, S is the stress matrix, and J is the Jacobian matrix. At this time, the internal force inside a single unit in the macro model is calculated.

[0112] Since then, in each process, the average stress of the entire micro model is solved in parallel using CPU multi-threading technology (0penMP) according to the stress of each Gaussian point on the micro model as the average stress on a single Gaussian point of a single unit of the macro model, and then the internal force of the macro model unit is calculated based on the average stress on the four Gaussian points of a single unit of the macro model.

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

[0114] S91, collecting internal force vectors;

[0115] After waiting for each process to complete the internal force calculation of all units in the process, all internal force vectors of each unit in each process are sent to process 0 through the MPI communication protocol. After waiting for the communication to be completed, that is, after process 0 completes the collection, the complete internal force vector of the macro model is obtained.

[0116] S92, calculate the macro model displacement;

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

[0118] KU=F

[0119] Since then, the internal force calculation of all units in each process has been completed. Process 0 collects the calculation results of the internal forces of all units in the macro model to obtain the internal force vector of the macro model, and calculates the displacement of the macro model based on the internal force vector.

[0120] The specific implementation method of step S10 is: set the external force vector F obtained by using the boundary condition in step S2 to a large number out The internal force vector F obtained with S9 in Add them together to get F r , use F r Solve the stiffness matrix to get the new displacement U n , which converges with the previous step.

[0121] If the calculation converges, F r The calculated U n This is the final displacement. If r does not meet the convergence condition, the displacement U obtained by the previous external force calculation is n As the displacement U of the new step n+1 , with U n+1 Set it as the boundary condition of the microscopic model in S6 and repeat the above steps until convergence.

[0122] The specific implementation method of step S11 is: after completing the calculation of the model displacement in S10, the stress and strain of each unit are calculated, and the calculation is performed using the method in S41, but the data is copied to the GPU for calculation, that is;

[0123] ε=BU

[0124] σ=Dε

[0125] Since this part only involves a large number of matrix operations, GPU parallelism, i.e. CUDA, is used for post-processing in the finite element method. First, the displacement matrix, B matrix and material matrix of all units are copied from the CPU memory to the GPU at one time, and the number of blocks (block) and threads (thread) are set in the GPU. For the calculation of strain, each thread calculates a row in the B matrix and a column in the U matrix, so that all threads can calculate the strain results of all units in parallel. For the calculation of stress, each thread calculates a row in the D matrix and a column in the strain matrix, so that all threads can calculate the stress results of all units in parallel. So far, the stress and strain at each Gaussian point of all units in the macro model have been obtained.

[0126] Since stress and strain are accurate solutions obtained at the integration points, and the strain and stress at the integration points are located inside the unit, in order to improve the global strain and stress results, the node values ​​and shape functions are used for linear interpolation to obtain the results at each node. The calculation formula is as follows:

[0127]

[0128] When multiple units share a node, the value at the node comes from the interpolation of several adjacent unit integration points, resulting in multiple stress and strain values ​​at the node, which requires smoothing. For example, for the strain and stress at a common node, the outer surface is a node shared by four units.

[0129] The calculation method is:

[0130]

[0131] Since then, we have completed the heterogeneous parallelization of CPU-GPU based on clusters and can effectively utilize computing resources to efficiently realize the simulation of anisotropic materials.

[0132] The beneficial effects of the present invention are as follows: the present invention realizes an efficient algorithm for various anisotropic materials based on the advantages of different parallel modes, firstly establishes a large-scale macro and micro structure and discretizes the grid, then uses the DDM technology to perform regional decomposition on the large-scale discrete model, and uses OpenMP to accelerate the loop to parallelly calculate the stiffness matrix of each unit. After iteratively solving the displacement of the macro model, the stress and strain on each Gaussian point of the macro model are calculated. The macro model unit is evenly divided into each process using the MPI technology, and the micro model boundary conditions are calculated using the macro model strain, and then the micro model displacement is calculated, and the average stress of the micro model is further calculated. The average stress of the micro model is used as the stress on a single Gaussian point of the macro model, and then the internal force of the macro model is calculated, and all the internal forces on each process are reduced using the MPI technology. The internal force is used to calculate the new displacement, and the residual is made with the original displacement and iterated, and the final displacement is obtained after the iteration converges. After the displacement calculation is completed, the CUDA of the GPU is used for post-processing parallelism. A large-scale heterogeneous parallel simulation calculation method based on regional decomposition and combined with multi-process, multi-threading and GPU parallelism is proposed by adopting modern large-scale computing. The idea of ​​"divide and conquer" is used to solve the problems of insufficient computing resources and long computing time, thereby realizing the calculation of large-scale anisotropic material structure response. BRIEF DESCRIPTION OF THE DRAWINGS

[0133] Figure 1 It is a flow chart of the method for solving large-scale cluster-CPU-GPU heterogeneous parallelism of the present invention;

[0134] Figure 2 Schematic diagram of macro and micro structures;

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

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

[0137] Figure 5This is a schematic diagram of K-way partitioning;

[0138] Figure 6 is the partition result diagram;

[0139] Figure 7 Collect and distribute schematic diagrams for MPI;

[0140] Figure 8 is the displacement result cloud map;

[0141] Fig. 9 It is a schematic diagram of shared nodes;

[0142] Fig.10 is the stress-strain cloud diagram; DETAILED DESCRIPTION

[0143] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the invention, not all of them. 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.

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

[0145] S1. Establish any large-scale macro-geometric model and micro-geometric model, discretize the macro-geometric model and micro-geometric model to be analyzed, divide the grid, and generate the required grid data;

[0146] The specific implementation steps are as follows: according to the complexity of macro-scale and micro-scale, a three-dimensional tetrahedral basic grid unit with reasonable density is established. The macro-model and micro-model are as follows: Figure 2 The model saves the number of model units and the nodes contained in each unit; it saves the number of nodes and the three-dimensional coordinates of all nodes.

[0147] S2. Set boundary conditions, constraints, and material parameters for the macro model, and use CPU multi-threading technology (OpenMP) to parallelly calculate the unit stiffness matrix of the macro model.

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

[0149] Set the material parameters of the macro model and micro model respectively. The macro model is regarded as an isotropic material, and the Poisson's ratio and elastic modulus are set according to the requirements. The micro model is an anisotropic material, and the elastic modulus and Poisson's ratio in each direction need to be set. For the macro model, it is also necessary to set boundary conditions and constraints to determine the location of the displacement boundary and fixed constraints.

[0150] S22, calculation of element stiffness matrix;

[0151] According to the set macro model material parameters, the macro model unit stiffness matrix is ​​calculated. First, the shape function of the three-dimensional tetrahedral unit is calculated. For the tetrahedral unit used, Figure 3 As shown, its volume coordinates are defined as;

[0152]

[0153] The shape function in the parametric coordinate system is equal to the corresponding volume coordinate. The four shape functions corresponding to the four nodes are:

[0154]

[0155] The partial derivative of the shape function in each direction of the local coordinate system is then multiplied by the node coordinates to obtain the Jacobian matrix. The method for solving the Jacobian matrix is ​​as follows;

[0156]

[0157] Multiply the inverse matrix of the Jacobi matrix by the partial derivative of the shape function to get the B matrix;

[0158]

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

[0160]

[0161] Among them, E and v are the elastic modulus and Poisson’s ratio of the macro model, respectively.

[0162] The expression of the unit stiffness matrix is ​​obtained by deducing the energy principle. Since the expression of the stiffness matrix is ​​an integral operation, the calculation of the unit stiffness matrix is ​​realized by introducing numerical integration. The results of hammer numerical integration are as follows;

[0163]

[0164] S23, OpenMP parallel computing;

[0165] Since the calculation of the element stiffness matrix of each element is independent, and the calculation of the element stiffness matrix is ​​only performed between a single loop, and involves fewer variables, the OpenMP multi-threading technology is used to perform parallel calculation of the element stiffness matrix. The elements of the macro model are divided into blocks, and each thread calculates the element stiffness matrix of the element assigned to each thread, such as Figure 4 shown.

[0166] S3. Use the regional decomposition algorithm to partition the model according to the set value, divide the grid model to be solved into multiple calculation sub-regions and make each partition carry all its own parameters. Iterate in each partition according to the partition information and solve the displacement field of the complete model. During the solution process, use the open source library STXXL to place the data in the memory on the hard disk, so as to reduce memory usage and quickly read the hard disk data.

[0167] S31, model domain decomposition;

[0168] The number of regions N to be divided is determined based on the calculated model size. The relationship between the number of regions to be divided and the partitioning method is usually that when N is less than 8, a multi-level recursive bisection method is used, and when N is greater than or equal to 8, a multi-level K-way partitioning method is used. The partitioning diagram is shown in the figure. Figure 5 shown.

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

[0170] S32, partition stiffness matrix calculation;

[0171] After completing the regional division of the model, the partition stiffness matrix calculation of each partition is integrated according to the partition situation. According to the partition unit information in S31, the unit stiffness matrix of each partition containing unit is extracted, and the partition stiffness matrix is ​​constructed according to the printed node degree of freedom encoding of the partition re-edited node information.

[0172] S33, iteratively solving the displacement field;

[0173] First, the global Schwarz solution iteration in the regional decomposition algorithm is:

[0174]

[0175] where ω is an arbitrary function that maps the domain Ω to R, and u nis an approximate solution to the Poisson problem. This formula is used to solve the subproblems in all subdomains of the regional algorithm at the same time, and its convergence is proved using the maximum principle. On this basis, the residual is introduced as follows;

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

[0177] The restricted additive Schwarz algorithm (RAS) on the continuous level can be defined.

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

[0179] KU=F

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

[0181]

[0182] Since the RAS algorithm is defined based on the continuous level, and 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, let D be the diagonal matrix of K, then;

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

[0184] Define the residual at the discrete level as r n =F-KU n , the Jacobi algorithm is divided into blocks, and the block Jacobi algorithm can be defined as;

[0185]

[0186] Among them, K ii is the overall stiffness matrix of each partition, is the displacement vector of the nth iteration step of the i-th partition. At this time, R is introduced i As N to N i The restriction operator, As N to N i The transpose operator of At this time, the block Jacobi algorithm can be written as follows:

[0187]

[0188] Based on the above formula, the RAS algorithm applicable to any number of partitions at the discrete level is:

[0189]

[0190] For the block Jacobi algorithm, the matrix expansion with the partition number of 2 is:

[0191]

[0192] in, Represents the displacement of the first partition at the nth iteration step The non-overlapping part of Represents the displacement of the first partition at the nth iteration step overlapping part. Represents the residual of the first partition in the nth iteration The non-overlapping part, Represents the overlapped residuals of the first and second partitions at the nth iteration.

[0193] The final linear algebraic form obtained by simplifying the above formula is:

[0194]

[0195] in, is the local correction of the discrete level solution in the nth step,

[0196] Each partition solves its own linear equations, and after obtaining the results, the overlapping parts are multiplied by one-half and added. This form is applicable to all partition numbers. When the number of partitions is greater than 2, the regional division situation can be rewritten with reference to the above formula. It is only necessary to reflect the partition situation of the overlapping nodes in the constant term matrix. At this point, the displacement field solution of the RAS algorithm based on regional decomposition is completed.

[0197] S34, memory data storage;

[0198] For the displacement field solution process mentioned above, for models with tens of millions of units, a large amount of information needs to be stored in the memory in actual calculations, such as unit number, node coordinates and other model information, as well as the B matrix, material matrix and other finite element model information of each unit. This information will occupy a large amount of memory in the calculation. The open source library STXXL can be used to store memory data in the hard disk, thereby reducing memory usage and quickly reading hard disk data.

[0199] This open source library uses the most popular parallel disk model (PDM). In this model, file I / O is explicitly handled by the program. The model divides the M-byte data in memory into D blocks and transfers them between D hard disks to amortize the delay. Usually, the interface of the external storage algorithm assumes that it reads input from the external storage, that is, the hard disk, and writes the output to the external storage. STXXL uses the idea of ​​pipeline to equip the external memory algorithm with a new interface, allowing them to provide the output as a data stream directly to the algorithm that uses the output instead of writing it to the memory first. In addition, STXXL uses synchronous I / O calls to run asynchronously in a separate thread, that is, for each file, there is a read request queue and a write request queue and a thread. The main thread sends the request to the file queue, and the thread associated with the file executes the request in FIFO order. Using this method, the data in the memory can be stored on the hard disk, thereby reducing the memory usage.

[0200] S4. According to the calculated displacement results, the strain and stress at each Gauss point of each unit of the macro model are calculated;

[0201] The specific implementation steps are: calculate the stress and strain at the Gaussian point according to the three-dimensional tetrahedron unit used. The three-dimensional tetrahedron unit has four Gaussian points. It is necessary to calculate the stress and strain in each direction at each Gaussian point respectively. The calculation method is as follows;

[0202] ε=BU

[0203] Among them, ε is the strain matrix at a single Gaussian point, recording the strains and shear strains in three directions, B is the B matrix at the Gaussian point, and U is the displacement matrix of each node on the unit sorted according to the unit connection order. After obtaining the strain matrix, the stress matrix is ​​calculated;

[0204] σ=Dε

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

[0206] Since then, the calculation of strain and stress at each Gaussian point of each unit in the macro model has been completed based on the calculated displacement results;

[0207] S5. Use the multi-process parallel technology (MPI) of the cluster to evenly distribute the units of the macro model 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 the displacement of the microscopic model at each Gaussian point needs to be calculated, and this involves more calculations rather than a single cycle, the MPI multi-process technology is used for parallel processing. First, the macro unit is divided into blocks, and the stress, strain and displacement of each Gaussian point on each unit are distributed to each process using the MPI communication protocol, such as Figure 7 shown.

[0209] S6. Reconstruct the model at each Gauss point of each unit of the macro model under each process, use the micro model as the microstructure of the macro model, and use the strain at the Gauss point to calculate the displacement boundary condition of the microstructure at the point;

[0210] The specific implementation steps are: reconstruct the model at each Gauss point on each unit of the macro model, and use the strain at the Gauss point to calculate the displacement boundary of the micro model. Since the finite element method first calculates the stress and strain at the Gauss point, and then uses the interpolation method to calculate the stress and strain at other positions of the unit, the micro model is regarded as a point on the macro model. Calculate the new displacement and stress and strain of the micro model at this time, and then regard the result as the new stress at the Gauss point to achieve the purpose of cross-scale calculation.

[0211] S7, using CPU multi-threading technology (OpenMP) in each process to parallelly calculate the unit stiffness matrix of each unit on the microstructure, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each unit of the micromodel according to the displacement results;

[0212] S71, microscopic model boundary calculation;

[0213] After distributing the macro model to each process according to the previous step, each process needs to perform a complete micro model calculation at each Gauss point in the unit. At this time, each process contains all the model information of the micro model and the strain at each Gauss point. First, extract the boundary nodes and coordinates of the micro model, and use the strain at the Gauss point multiplied by the coordinate to obtain the displacement boundary conditions of the micro model;

[0214] Δd=εd

[0215] Where Δd is the displacement of each boundary node on the microscopic model, and d is the coordinate of each boundary node on the microscopic model. At this time, the boundary conditions of the microscopic model are obtained.

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

[0217] The same finite element theory in S22 is used to calculate the shape function and other parameters in the micro-model unit stiffness matrix. Since the macro-model is regarded as an isotropic material and the micro-model is regarded as an anisotropic material, the calculation method of the micro-model is different from that of the macro-model when calculating the material matrix D. At this time, the calculation method of the material matrix D of the micro-model is:

[0218]

[0219] According to the different forms of anisotropic materials, 36 different parameters in the matrix are determined. And the stiffness matrix of each element is calculated in the same way as in S2.

[0220] S73, OpenMP parallel computing;

[0221] According to the same method in S23, OpenMP is used to accelerate the calculation of the element stiffness matrix and realize the simultaneous scheduling of processes and threads.

[0222] S74, overall stiffness matrix assembly;

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

[0224] S75, solution of displacement, stress and strain;

[0225] According to the same method in S41, the displacement field in the microscopic model and the stress-strain value at each Gaussian point on the microscopic model are solved at this time.

[0226] S8. In each process, the average stress of the entire micro-model is solved in parallel using CPU multi-threading technology (OpenMP) according to the stress of each Gaussian point on the micro-model as the average stress of a single Gaussian point of a single unit of the macro-model, and then the internal force of the unit of the macro-model is calculated according to the average stress of four Gaussian points of a single unit of the macro-model;

[0227] S81, calculation of average stress of micro model;

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

[0229]

[0230] in, is the average stress of the model, σ ijis the stress at each Gauss point.

[0231] S82, OpenMP parallelism;

[0232] Since the average stress calculation of each unit in the microscopic model is independent of each other, the S23 method is used to accelerate the calculation of the average stress, and OpenMP is used to enable multi-threading to accelerate the loop.

[0233] S83, calculation of internal forces of micro-model units;

[0234] After obtaining the average stress of a single microscopic model, the macroscopic model is returned, that is, the stress at a single Gaussian point of the macroscopic model is equivalent to the average stress of the microscopic model at the Gaussian point, and the internal force of the macroscopic model at this time is calculated based on the average stress;

[0235]

[0236] Among them, B is the B matrix of the macro unit, S is the stress matrix, and J is the Jacobian matrix. At this time, the internal force inside a single unit in the macro model is calculated.

[0237] S9, each process completes the internal force calculation of all units in the process, process 0 collects the calculation results of the internal forces of all units of the macro model to obtain the internal force vector of the macro model, and calculates the displacement of the macro model according to the internal force vector;

[0238] S91, collecting internal force vectors;

[0239] After waiting for each process to complete the internal force calculation of all units in the process, all internal force vectors of each unit in each process are sent to process 0 through the MPI communication protocol. After waiting for the communication to be completed, that is, after process 0 completes the collection, the complete internal force vector of the macro model is obtained.

[0240] S92, calculate the macro model displacement;

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

[0242] KU=F

[0243] S10, iterating the displacement judgment residual obtained by using the displacement calculated by the displacement boundary of the macro model and the internal force vector calculation, and obtaining the displacement result of the macro model after judging convergence;

[0244] The specific implementation steps are: set the external force vector F obtained by using the boundary conditions in S2 to a large number out The internal force vector F obtained with S9 in Add them together to get F r , use F rSolve the stiffness matrix to get the new displacement U n , and the previous step determines convergence. If the calculation converges, F r The calculated U n This is the final displacement. If r does not meet the convergence condition, the displacement U obtained by the previous external force calculation is n As the displacement U of the new step n+1 , with U n+1 Set it as the boundary condition of the microscopic model in S6, repeat the above steps until convergence, and the displacement result is as follows Figure 8 shown.

[0245] S11. Calculate the stress and strain of the macro model based on the displacement results through GPU, quickly implement post-processing of the calculation, and realize the cross-scale heterogeneous parallel algorithm of the anisotropic material composition structure.

[0246] The specific implementation steps are:

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

[0248] The specific implementation steps are: after completing the calculation of the model displacement in S10, the stress and strain of each unit are calculated, and the calculation is performed using the method in S41, but the data is copied to the GPU for calculation, that is;

[0249] ε=BU

[0250] σ=Dε

[0251] Since this part only involves a large number of matrix operations, GPU parallelism, i.e. CUDA, is used for post-processing in the finite element method. First, the displacement matrix, B matrix and material matrix of all units are copied from the CPU memory to the GPU at one time, and the number of blocks (block) and threads (thread) are set in the GPU. For the calculation of strain, each thread calculates a row in the B matrix and a column in the U matrix, so that all threads can calculate the strain results of all units in parallel. For the calculation of stress, each thread calculates a row in the D matrix and a column in the strain matrix, so that all threads can calculate the stress results of all units in parallel. So far, the stress and strain at each Gaussian point of all units in the macro model have been obtained.

[0252] Since stress and strain are accurate solutions obtained at the integration points, and the strain and stress at the integration points are located inside the unit, in order to improve the global strain and stress results, the node values ​​and shape functions are used for linear interpolation to obtain the results at each node. The calculation formula is as follows:

[0253]

[0254] When multiple units share a node, the value at the node comes from the interpolation of several adjacent unit integration points, resulting in multiple stress and strain values ​​at the node, which requires smoothing. For the strain and stress at the common node, the outer surface is a node shared by four units, such as Fig. 9 As shown:

[0255] The calculation method is:

[0256]

[0257] Since then, the heterogeneous parallelization of CPU-GPU based on cluster has been completed, and the simulation of anisotropic materials can be efficiently realized by effectively utilizing computing resources. The post-processing results are as follows: Fig.10 shown.

[0258] The present invention applies the DDM method and heterogeneous parallel computing strategy to the finite element calculation of anisotropic materials. According to the advantages of different parallel methods, the numerical computing ability of the computer is fully utilized to realize the calculation of large-scale models, and the rapid calculation of large-scale anisotropic material structure response can be realized.

[0259] The above description does not limit the present invention in any form. Although the present invention has been disclosed by the above embodiments, it is not used to limit the present invention. Any technician familiar with the profession can make changes or modifications to equivalent embodiments of equivalent changes by using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent changes and modifications made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural responses, characterized in that: The following steps are involved: S1. Establish any large-scale macro-geometric model and micro-geometric model, discretize the macro-geometric model and micro-geometric model to be analyzed, divide the grid, and generate the required grid data; S2, setting boundary conditions, constraints, and material parameters for the macro model, and using CPU multi-threading technology (OpenMP) to parallelly calculate the unit stiffness matrix of the macro model; S3, using the regional decomposition algorithm to partition the model according to the set value, dividing the solution grid model into multiple calculation sub-regions while making each partition carry all its own parameters; Iterate in each partition according to the partition information and solve the displacement field of the complete model. During the solution process, use the open source library STXXL to place the data in the memory on the hard disk to reduce memory usage and quickly read hard disk data. S4. According to the calculated displacement results, the strain and stress at each Gauss point of each unit of the macro model are calculated; S5. Use the multi-process parallel technology (MPI) of the cluster to evenly distribute the units of the macro model 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 Gauss point of each unit of the macro model under each process, use the micro model as the microstructure of the macro model, and use the strain at the Gauss point to calculate the displacement boundary condition of the microstructure at the point; S7, using CPU multi-threading technology (OpenMP) in each process to parallelly calculate the unit stiffness matrix of each unit on the microstructure, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each unit of the micromodel according to 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) according to the stress of each Gaussian point on the micro-model as the average stress of a single Gaussian point of a single unit of the macro-model, and then the internal force of the unit of the macro-model is calculated according to the average stress of four Gaussian points of a single unit of the macro-model; S9, each process completes the internal force calculation of all units in the process, process 0 collects the calculation results of the internal forces of all units of the macro model to obtain the internal force vector of the macro model, and calculates the displacement of the macro model according to the internal force vector; S10, iterating the displacement judgment residual obtained by using the displacement calculated by the displacement boundary of the macro model and the internal force vector calculation, and obtaining the displacement result of the macro model after judging convergence; S11. Calculate the stress and strain of the macro model based on the displacement results through GPU, quickly implement post-processing of the calculation, and realize the cross-scale heterogeneous parallel algorithm of the anisotropic material composition structure.

2. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S1 is: establishing a three-dimensional tetrahedral basic grid unit with a reasonable density according to the complexity of the macro scale and the micro scale; the model saves the number of model units and the nodes contained in each unit; and saves the number of nodes and the three-dimensional coordinates of all nodes.

3. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S2 is: S21, setting material parameters, boundary conditions and constraints; Set the material parameters of the macro model and micro model respectively. The macro model is regarded as an isotropic material, and the Poisson's ratio and elastic modulus are set according to the requirements. The micro model is an anisotropic material, and the elastic modulus and Poisson's ratio in each direction need to be set. For the macro model, boundary conditions and constraints need to be set to determine the positions of the displacement boundaries and fixed constraints. S22, calculation of element stiffness matrix; According to the set macro model material parameters, the macro model unit stiffness matrix is ​​calculated. First, the shape function of the three-dimensional tetrahedral unit is calculated. For the tetrahedral unit used, its volume coordinate is defined as; The shape function in the parametric coordinate system is equal to the corresponding volume coordinate. The four shape functions corresponding to the four nodes are: The partial derivative of the shape function in each direction of the local coordinate system is then multiplied by the node coordinates to obtain the Jacobian matrix. The method for solving the Jacobian matrix is ​​as follows; Multiply the inverse matrix of the Jacobi matrix by the partial derivative of the shape function to get the B matrix; The D matrix of an isotropic three-dimensional tetrahedral element is as follows; in, E and v are the elastic modulus and Poisson’s ratio of the macro model, respectively; The expression of the unit stiffness matrix is ​​obtained by deducing the energy principle. Since the expression of the stiffness matrix is ​​an integral operation, the calculation of the unit stiffness matrix is ​​realized by introducing numerical integration. The results of hammer numerical integration are as follows; S23, OpenMP parallel computing; Since the calculation of the element stiffness matrix of each element is independent, and the calculation of the element stiffness matrix is ​​only performed between a single loop, and involves fewer variables, the OpenMP multi-threading technology is used to perform parallel calculation of the element stiffness matrix; the elements of the macro model are divided into blocks, and each thread calculates the element stiffness matrix of the element assigned to each thread; At this point, the setting of boundary conditions, constraints, and material parameters for the macro model has been completed, and the unit stiffness matrix of the macro model is calculated in parallel using CPU multi-threading technology (0penMP).

4. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S3 is: S31, model domain decomposition; The number of divided regions N is determined according to the calculated model size. The relationship between the number of divided regions and the partitioning method is usually that when N is less than 8, the multi-level recursive bisection method is used, and when N is greater than or equal to 8, the multi-level K-way partitioning method is used; After the partitioning is completed, the nodes and units of each partition are renumbered, and the node information and unit information of each partition before and after the renumbering as well as the unit information at the boundary of each partition are saved; At this time, the original finite element model is divided into N blocks of regions Ω0~Ω N-1 S32, partition stiffness matrix calculation; After completing the regional division of the model, the partition stiffness matrix calculation of each partition is integrated according to the partition situation; according to the partition unit information in S31, the unit stiffness matrix of each partition containing unit is extracted, and the partition stiffness matrix is ​​constructed according to the node freedom degree coding printed according to the partition re-edited node information; S33, iteratively solving the displacement field; First, the global Schwarz solution iteration in the regional decomposition algorithm is: where ω is an arbitrary function that maps the domain Ω to R, and u n is an approximate solution to the Poisson problem; this formula is used to simultaneously solve the subproblems in all subdomains of the regional algorithm, and its convergence is proved using the maximum principle; on this basis, the residual is introduced as follows; r n =f+Δ(u n ) That is, the restricted additive Schwarz algorithm (RAS) on the continuous level can be defined; In the finite element elastic system, for quasi-static solution, according to Hooke's law and the principle of minimum potential energy: KU=F Among them, K is the global stiffness matrix, F is the force matrix, and U is the displacement vector result; based on the regional decomposition algorithm, the above linear system can be converted into the following block matrix form: Since the RAS algorithm is defined based on the continuous level, and 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; introduce the block Jacobi algorithm, let D be the diagonal matrix of K, then; YOU n+1 =YOU n +(F-KU n ) Define the residual at the discrete level as r n =F-KU n , the Jacobi algorithm is divided into blocks, and the block Jacobi algorithm can be defined as; Among them, K ii is the overall stiffness matrix of each partition, is the displacement vector of the nth iteration step of the i-th partition; at this time, R is introduced i As N to N i The restriction operator, As N to N i The transpose operator of At this time, the block Jacobi algorithm can be written as follows: Based on the above formula, the RAS algorithm applicable to any number of partitions at the discrete level is: For the block Jacobi algorithm, the matrix expansion with the partition number of 2 is: in, Represents the displacement of the first partition at the nth iteration step The non-overlapping part of Represents the displacement of the first partition at the nth iteration step The overlapping part of Represents the residual r1 of the first partition in the nth iteration n The non-overlapping part, Represents the overlapped residuals of the first partition and the second partition at the nth iteration step; The final linear algebraic form obtained by simplifying the above formula is: in, is the local correction of the discrete level solution in the nth step, Each partition solves its own linear equations, and after obtaining the results, the overlapping parts are multiplied by one-half and added; this form is applicable to all partition numbers. When the number of partitions is greater than 2, the regional division situation can be rewritten with reference to the above formula. It is only necessary to reflect the partition situation of the overlapping nodes in the constant term matrix; thus, the displacement field solution of the RAS algorithm based on regional decomposition is completed; S34, memory data storage; For the above displacement field solution process, for models with tens of millions of units, a large amount of information needs to be stored in the memory in actual calculations, such as unit number, node coordinates and other model information, as well as the B matrix, material matrix and other finite element model information of each unit; this information will occupy a large amount of memory in the calculation, and the open source library STXXL can be used to store memory data in the hard disk, thereby reducing memory usage and quickly reading hard disk data; This open source library uses the most popular parallel disk model (PDM); in this model, file I / O is explicitly handled by the program. The model divides the M-byte data in the memory into D blocks and transmits them between D hard disks to spread the delay; usually the interface of the external storage algorithm assumes that it reads input from the external storage, that is, the hard disk, and writes the output to the external storage; STXXL uses the idea of ​​pipeline to equip the external memory algorithm with a new interface, allowing them to provide the output as a data stream directly to the algorithm that uses the output instead of writing it to the memory first; in addition, STXXL uses synchronous I / O calls to run asynchronously in a separate thread, that is, for each file, there is a read request queue and a write request queue and a thread; the main thread sends the request to the file queue, and the thread associated with the file executes the request in FIFO order; this method can be used to store data in memory on the hard disk to reduce memory usage; At this point, we have completed the use of the regional decomposition algorithm to partition the model according to the set values, dividing the grid model to be solved into multiple calculation sub-regions while allowing each partition to carry all of its own parameters; according to the partition information, iterate in each partition and solve the displacement field of the complete model. During the solution process, the open source library STXXL is used to place the data in the memory on the hard disk to reduce memory usage and quickly read hard disk data.

5. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S4 is: S41. Solve the stress and strain at the Gauss point; The stress and strain at the Gaussian point are calculated based on the three-dimensional tetrahedral unit used. The three-dimensional tetrahedral unit has four Gaussian points. 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 strains in three directions and the shear strains in three directions, B is the B matrix at the Gaussian point, and U is the displacement matrix of each node on the unit sorted in the order of unit connection; after obtaining the strain matrix, the stress matrix is ​​calculated; σ=Dε Where σ is the stress matrix at a single Gaussian point, recording the stresses and shear stresses in three directions, and the D matrix is ​​the material matrix of the unit; Since then, the calculation of strain and stress at each Gaussian point of each unit in the macro model has been completed based on the calculated displacement results.

6. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S5 is: since it is necessary to calculate the displacement of the micro model at each Gauss point, and more calculations are involved rather than a single cycle, MPI multi-process technology is used for parallelization; first, the macro unit is divided into blocks, and the stress, strain and displacement of each Gauss point on each part of the unit are distributed to each process using the MPI communication protocol.

7. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S6 is: reconstruct the model at each Gauss point on each unit of the macro model, and use the strain at the Gauss point to calculate the displacement boundary of the micro model; since the finite element method first calculates the stress and strain at the Gauss point, and then uses the interpolation method to calculate the stress and strain at other positions of the unit, the micro model is regarded as a point on the macro model; calculate the new displacement and stress and strain of the micro model at this time and then regard the result as the new stress at the Gauss point to achieve the purpose of cross-scale calculation.

8. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S7 is: S71, microscopic model boundary calculation; After distributing the macro model to each process according to the previous step, each process needs to perform a complete micro model calculation at each Gauss point in the unit; at this time, each process contains all the model information of the micro model and the strain at each Gauss point; first extract the boundary nodes and coordinates of the micro model, and use the strain at the Gauss point multiplied by the coordinate to obtain the displacement boundary conditions of the micro model; Δd=εd Where Δd is the displacement of each boundary node on the microscopic model, and d is the coordinate of each boundary node on the microscopic model; at this time, the boundary conditions of the microscopic model are obtained; S72, micro-model element stiffness matrix calculation; The same finite element theory in S22 is used to calculate the shape function of the micro-model unit stiffness matrix. Since the macro-model is regarded as an isotropic material and the micro-model is regarded as an anisotropic material, the calculation method of the micro-model is different from that of the macro-model when calculating the material matrix D. At this time, the calculation method of the material matrix D of the micro-model is: Determine 36 different parameters in the matrix according to different forms of anisotropic materials; and calculate the stiffness matrix of each element in the same way as in S2; S73, OpenMP parallel computing; According to the same method as in S23, OpenMP is used to accelerate the calculation of the element stiffness matrix and realize the simultaneous scheduling of processes and threads; S74, overall stiffness matrix assembly; After the OpenMP calculation of the unit stiffness matrix is ​​completed, the unit stiffness matrix is ​​assembled into the global stiffness matrix in sequence according to the degree of freedom index of each node in each unit to complete the assembly of the unit stiffness matrix and obtain the global stiffness matrix K; S75, solution of displacement, stress and strain; According to the same method as in S41, the displacement field in the micro-model and the stress-strain values ​​at each Gaussian point on the micro-model are solved at this time; At this point, the displacement, stress and strain of the micro model at a single Gauss point have been solved, and the CPU multi-threading technology (0penMP) is used in each process to parallelly calculate the unit stiffness matrix of each unit in the micro structure, assemble the overall stiffness matrix and solve the displacement field. According to the displacement results, the strain and stress at each Gauss point of each unit in the micro model are calculated.

9. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S8 is: S81, calculation of average stress of micro model; According to the stress on all Gaussian points on a single microscopic model calculated in S75, the average stress on the microscopic model at the Gaussian point is calculated. Based on the homogenization theory of simulating macroscopic mechanical properties of the composite material mesoscopic unit cell model, the average stress is calculated as follows: in, is the average stress of the model, σ ij is the stress at each Gauss point; S82, OpenMP parallelism; Since the average stress calculation of each unit in the microscopic model is independent of each other, the S23 method is used to accelerate the calculation of the average stress, and OpenMP is used to enable multi-threading to accelerate the loop; S83, calculation of internal forces of micro-model units; After obtaining the average stress of a single microscopic model, the macroscopic model is returned, that is, the stress at a single Gaussian point of the macroscopic model is equivalent to the average stress of the microscopic model at the Gaussian point, and the internal force of the macroscopic model at this time is calculated based on the average stress; Among them, B is the B matrix of the macro unit, S is the stress matrix, and J is the Jacobi matrix; at this time, the internal force inside a single unit in the macro model is calculated; Since then, in each process, the average stress of the entire micro model is solved in parallel using CPU multi-threading technology (0penMP) according to the stress of each Gaussian point on the micro model as the average stress on a single Gaussian point of a single unit of the macro model, and then the internal force of the macro model unit is calculated based on the average stress on the four Gaussian points of a single unit of the macro model.

10. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S9 is: S91, collecting internal force vectors; After waiting for each process to complete the internal force calculation of all units in the process, all internal force vectors of each unit in each process are sent to process 0 through the MPI communication protocol; after waiting for the communication to be completed, that is, after process 0 completes the collection, the complete internal force vector of the macro model is obtained; S92, calculate the macro model displacement; Get the internal force vector F in Then, the new displacement field U is calculated using the following formula: new ; KU=F Since then, the internal force calculation of all units in each process has been completed. Process 0 collects the calculation results of the internal forces of all units in the macro model to obtain the internal force vector of the macro model, and calculates the displacement of the macro model based on the internal force vector.

11. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S10 is: the external force vector F obtained by using the boundary conditions in S2 is set to a large number out The internal force vector F obtained by S9 in Add them together to get F r , use F r Solve the stiffness matrix to get the new displacement U n , converges with the previous step; If the calculation converges, F r The calculated U n is the final displacement; if r does not meet the convergence condition, the displacement U obtained by the previous external force calculation is n As the displacement U of the new step n+1 , with U n+1 Set it as the boundary condition of the microscopic model in S6 and repeat the above steps until convergence.

12. According to the cluster-CPU-GPU heterogeneous parallel method for solving large-scale anisotropic structural response according to claim 1, the specific implementation method of step S11 is: after completing the calculation of the model displacement in S10, the stress and strain of each unit are calculated, and the calculation is performed using the method in S41, but the data is copied to the GPU for calculation, that is; ε=BU σ=Dε Since this part only involves a large number of matrix operations, GPU parallelization, i.e. CUDA, is used for post-processing in the finite element method. First, the displacement matrix, B matrix and material matrix of all units are copied from the CPU memory to the GPU at one time, and the number of blocks and threads is set in the GPU. For strain calculation, each thread calculates a row in the B matrix and a column in the U matrix, so that all threads can calculate the strain results of all units in parallel. For stress calculation, each thread calculates a row in the D matrix and a column in the strain matrix, so that all threads can calculate the stress results of all units in parallel. So far, the stress and strain at each Gaussian point of all units in the macro model are obtained. 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 unit, in order to improve the global strain and stress results, linear interpolation is performed using the node values ​​and shape functions to obtain the results at each node; the calculation formula is as follows: When multiple units share a node, the value at the node comes from the interpolation of several adjacent unit integration points, resulting in multiple stress and strain values ​​at the node, which requires smoothing. For example, for the strain and stress at a common node, the outer surface is a node shared by four units. The calculation method is: Since then, we have completed the heterogeneous parallelization of CPU-GPU based on clusters and can effectively utilize computing resources to efficiently realize the simulation of anisotropic materials.

Images