CN106529011B - A kind of Parallel districts implementation method for SPH algorithm - Google Patents
A kind of Parallel districts implementation method for SPH algorithm Download PDFInfo
- Publication number
- CN106529011B CN106529011B CN201610968969.7A CN201610968969A CN106529011B CN 106529011 B CN106529011 B CN 106529011B CN 201610968969 A CN201610968969 A CN 201610968969A CN 106529011 B CN106529011 B CN 106529011B
- Authority
- CN
- China
- Prior art keywords
- particle
- parallel districts
- boundary
- parallel
- coordinate
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/30—Circuit design
- G06F30/36—Circuit design at the analogue level
- G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
Landscapes
- Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Microelectronics & Electronic Packaging (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a kind of Parallel districts implementation methods for SPH algorithm comprising following steps: creating the simulation model of pending simulation analysis and carries out Initialize installation；The pressure value, stress tensor and pseudo-viscosity value of each particle corresponding to phantom；The simulation model is scanned for analyzing in each set time step；The searching analysis includes each Parallel districts corresponding to determining simulation model and each Parallel districts boundary；State equation corresponding to phantom and governing equation；Simulation result is exported based on set termination condition.Present invention reduces memory spaces required for simulation process, and greatly reduce the traffic, and then greatly improve computational efficiency, and effectively reduce calculating cost.
Description
Technical field
The present invention relates to computer simulation technique field, particularly relate to realize for the Parallel districts of SPH algorithm
Method.
Background technique
Usually our object of which movement for being concerned about and the rule of deformation can be described using the differential equation, claim this differential side
Journey is called physical control equation, actually most moving object be all it is complex-shaped, physical control equation must be by
Numerical value calculating is carried out in computer help, is also numerical simulation or numerical simulation, is meant according to objective law using numerical value
The method simulation object movement under certain condition of calculating and deformation, this belongs to a kind of prediction, under present technical conditions
Prediction result and actual conditions have been able to closely.Thus, numerical simulation calculates in scientific research and engineer application
Have become the third research method for being listed in laboratory facilities and theoretical means, with visual result, repeatability is high, item
Part requires the advantages such as low and at low cost.
In general, numerical computation method can be divided into two kinds: another then be one is the numerical computation method based on grid
Mesh free numerical computation method；Wherein the numerical computation method based on grid is divided into Lagrangian method and Euler's method: institute
Stating Lagrangian method is that grid is allowed to change with the movement and deformation of object, and this method is used primarily in the numerical value meter of solid
It calculates, can accurately capture the change in shape of solid, irregular contour can be distributed to describe with irregular grid, be counted
It is high-efficient, it is adaptable, it is widely used in Solid Mechanics calculating, also there is a large amount of application in machinery industry；But it is positive because
Change for the grid of Lagrangian method with deformation, when deforming very big, distortion of the mesh is very serious, at this time cannot
Meet the basic demand for having grid values calculation method, therefore will cause time step reduction, the time-consuming increase calculated each time,
Even fail because mesh distortion seriously causes to calculate, limits this method application；The Euler's method is that grid is allowed to fix
Constant, the result obtained each time is all the physical quantity of the fixed position in space, is mainly used for the calculating of flow problem, because of grid
Immobilize, therefore will not influence to calculate because of distortion of the mesh, exactly because but fixed grid, therefore its cannot describe it is solid
The deformation of body, and do not have to confirm the measure at interface between multiple fluid in for fluid calculation, it is therefore desirable to it is additional to propose
Confirm the method for fluid boundary, Euler's method effectively handles the side at interface naturally still not as Lagrangian method at present
Method, therefore the method at its processing interface used in a fluid can only meet the requirement of calculating to a certain extent, not very
Accurately, it can satisfy fluid calculation requirement substantially, but allowed for solid calculating and be not able to satisfy accurately retouching for solid boundaries
It states.
From the point of view of two kinds have the characteristics of grid method, there is grid method all to have the defects that some not applicable, although can be with
It is complementary to one another, but cannot be used together due to the difference of description method, if the numerical value calculating that meet in Solid Mechanics is wanted
It asks, it is necessary to remove grid, meet the requirement of hydrodynamics method, it is necessary to which a kind of Lagrangian method describes more substance fluids
Interface.
And the SPH method based on non-mesh method just can satisfy these requirements, SPH method full name is smooth particle flux
Body dynamics method is a kind of non-mesh method of Lagrangian format, both can be used for the calculating of Solid Mechanics numerical value, can also
To be calculated for hydrodynamics numerical value；It is when being used for Solid Mechanics, because without grid, the problem of mesh distortion is not deposited
It is SPH method, when being used for hydrodynamics, because being Lagrangian method, does not need additional interface tracking technology
The interface of multiple fluid can be described accurately.
With going deep into for research, SPH method is gradually applied to underwater explosion and high velocity impact, in ocean wave simulation and
Also there is successful application in the simulation of slosh, be widely used in the research of spacecraft space fragment protective in recent years.But with
SPH method be gradually used in actual engineer application, existing SPH method also suffers from certain drawbacks, and specifically includes
Following problems:
For having for grid method, the relationship between node can be determined by grid, and be kept not in calculating process
Become, the grid of each subregion and node relationships calculate start before it has been determined that therefore parallel processing it is very not difficult, but
Be that correlation for non-mesh method, between particle is all changing all the time, thus each time before search all
It needs to repartition or adjust parallel regions, therefore efficiently judges the direction of subregion and the size and location of each subregion
Just become critically important.
Current existing SPH parallel method includes the parallel method of shared storage and the parallel method of distribution storage, wherein
The advantages of shared storage needs all cpu that can share same memory, this method is to program very simple, does not need to change
Code, the disadvantage is that scalability is very poor, it is high to hardware requirement；The parallel method of distributed storage is then to utilize network connection not
Same host, each host is independently to be calculated, and carries out necessary communication by network to realize that simulation is same
The shortcomings that different piece of problem, distributed storage is that program becomes very complicated, and advantage is that scalability is fine, therefore distribution is simultaneously
Capable program can be run on the computer of shared storage, but the concurrent program of shared storage cannot be used for distributed computing
Machine, therefore the SPH program for developing distributed parallel has better application prospect.Current distributed parallel SPH method is being searched
The parallel method of use is all simplified during rope, and the most commonly used is all give all particle informations to each to calculate section
Point, only each calculate node, which is responsible for searching for, oneself is responsible for part particle, therefore can give each node largely extra letter
Breath, it is thus possible to bring a large amount of storage burden, while calculate the whole for terminating that each calculate node is required to be responsible for every time
Particle information, which communicates, gives each calculate node, causes efficiency extremely low.It is searched parallel seeing as existing SPH in summary
The technical restriction of Suo Fangfa, the prior art can accomplish to give particle into each calculate node respectively, but cannot be guaranteed each
Relationship is searched between the load balancing of node and clearly determining node.
Summary of the invention
In view of defects in the prior art, the invention aims to provide a kind of Parallel districts for SPH algorithm
Implementation method to save memory space required for simulation process, and greatly reduces the traffic, and then greatly improves calculating effect
Rate.
To achieve the goals above, technical solution of the present invention:
A kind of Parallel districts implementation method for SPH algorithm, which comprises the steps of:
Step 1, the simulation model of the pending simulation analysis of creation simultaneously carry out Initialize installation to the simulation model created；
Step 2, the pressure value for calculating each particle corresponding to the simulation model after initialized setting, stress tensor with
And pseudo-viscosity value；
Step 3 scans for analyzing in each set time step to the simulation model, by determination it is each in terms of
Calculate neighbour's particle corresponding to particle；The searching analysis includes determining and storing each Parallel districts corresponding to simulation model
And each Parallel districts boundary；
State equation corresponding to step 4, phantom and governing equation；
Step 5 is based on set termination condition, exports corresponding simulation result.
Further, as of the invention preferred
Step 3 includes that step 31 determines each Parallel districts and each Parallel districts boundary corresponding to simulation model, packet
It includes:
Step 311 first recycles all particles corresponding to simulation model, is found respectively in three-dimensional space i.e. X
Direction, Y-direction, in Z-direction, coordinate smallest particles I and coordinate maximum particle II；Secondly by the coordinate and grain of above-mentioned particle II
The coordinate of sub- I subtracts each other, and obtains the maximum span value on tri- directions X, Y, Z；Finally by will be corresponding to the simulation model
The coordinate of each particle successively subtracts each other with the coordinate of particle I and divided by the maximum span value on the coordinate direction obtained, really
Determine and store the opposite dimensionless coordinate (X of each particlej-Xi/X0,Yj-Yi/Y0,Zj-Zi/Z0), wherein Xj、Yj、ZjTable respectively
Show X-direction, Y-direction, the particle coordinate in Z-direction；Xi、Yi、ZiRespectively indicate that X-direction, Y-direction, coordinate is the smallest in Z-direction
The coordinate of particle I, X0Indicate the maximum span value in X-direction, Y0Indicate the maximum span value in Y-direction, Z0It indicates in Z-direction
Maximum span value, j >=1；
Step 312 calculates separately square of the dimensionless coordinate in X-direction, Y-direction, Z-direction corresponding to all particles
With, and the vector direction that the smallest dimensionless coordinate of quadratic sum is constituted is determined as principal direction；
Step 313, based on identified principal direction, find out the relative coordinate of each particle in a main direction, then according to
Relative coordinate is ranked up all particles from as low as big sequence；
Step 314, determined based on Parallel districts quantity to be delimited primary number corresponding to each Parallel districts with
And boundary particle serial number corresponding to each Parallel districts boundary position；
Step 315 determines each Parallel districts boundary position and stores the corresponding boundary of identified boundary particle serial number
Coordinate；
Step 316 delimited Parallel districts and be ranked up to the particle inside each Parallel districts, and pass through each parallel point
Area communicates to determine the number of particles situation of change in each Parallel districts, and based on identified number of particles situation of change tune
Whole Parallel districts boundary position, to reach load balancing.
Further, as of the invention preferred
All particles are ranked up using the bucket sort method of floating number in the step 313.
Further, as of the invention preferred
The step 314 includes:
Step 3141 obtains Parallel districts quantity corresponding to simulation model；
Step 3142 is divided by with identified Parallel districts quantity according to all populations corresponding to simulation model and is obtained
The integer obtained adds one to determine population initial inside each Parallel districts；
Step 3143 determines boundary particle serial number corresponding to each Parallel districts boundary position, i.e., parallel according to first
Population possessed by first Parallel districts of boundary particle serial number corresponding to subregion and the second Parallel districts boundary position,
And the n-th -2 Parallel districts of boundary particle serial number and (n-1)th Parallel districts in (n-1)th parallel area and n-th of parallel area
Boundary particle serial number plus (n-1)th parallel area population strategy, all Parallel districts are recycled, with determination
Boundary particle serial number corresponding to each Parallel districts boundary position, wherein n >=3.
Further, as of the invention preferred
The step 315 includes:
It is traversed according to the sequence that step 313 sequences, when traversing the boundary particle serial number of some Parallel districts,
The number of particle corresponding to the boundary particle serial number is extracted, and finds the opposite seat of the corresponding particle of the number in a main direction
Mark, this relative coordinate is the boundary position of present parallel subregion Yu next Parallel districts；Complete time of all particle sequences
It goes through to obtain the coboundary of all Parallel districts and lower boundary, wherein coboundary is that current Parallel districts are divided parallel with previous
The boundary coordinate in area, lower boundary are the boundary coordinate of itself and latter Parallel districts；And first Parallel districts coboundary is main
Particle coordinate corresponding to the initial position of direction, the last one Parallel districts lower boundary are grains corresponding to principal direction final position
Subcoordinate.
Further, as of the invention preferred
The step 316 includes:
Step 3161 delimited Parallel districts, and had from as low as big strategy to each Parallel districts according to principal direction
Particle be ranked up, with particle sequence corresponding to the current Parallel districts of determination；
Step 3162 determines that the number of particles in present parallel subregion changes feelings after time step iterative cycles each time
Condition: i.e. since the first Parallel districts, judge whether particle corresponding in current region is more than primary number, be then from grain
The particle having more is given into next Parallel districts in subsequence, and updates coboundary and the present parallel subregion of next Parallel districts
Lower boundary；
Step 3163 repeats step 3162, to be sequentially completed number of particles situation of change corresponding to each Parallel districts,
And up-and-down boundary position in Parallel districts is adjusted based on identified number of particles situation of change.
Compared with prior art, beneficial effects of the present invention:
Firstly, the present invention is efficiently solved existing by the scheme based on identified principal direction setting Parallel districts
Concurrent technique cannot be guaranteed the defect that relationship is searched between the load balancing of each calculate node and specific node, so that each meter
Particle is suitable on operator node, and the complexity searched for is O (N), this just can substantially ensure that and search to each calculate node
Quite, to realize load balancing substantially, and each calculate node is only needed to save and is needed the time-consuming that rope and governing equation calculate
Particle to be processed, it is not necessary to save the particle of other calculate nodes processing；The grain for enabling to each calculate node responsible simultaneously
Have the communication between clearly simple boundary and subregion very clear between son, to realize saving memory space, greatly
The traffic is reduced, the purpose of design of computational efficiency is greatly improved.
Detailed description of the invention
Fig. 1 is the corresponding process step figure of the Parallel districts implementation method of the present invention for SPH algorithm；
Fig. 2 is the corresponding process step example diagram of the Parallel districts implementation method of the present invention for SPH algorithm；
Fig. 3 is the example diagram that computational domain of the present invention is divided into 4 pieces of impartial regions；
Fig. 4 is the example diagram that the size of region of search of the present invention is greater than current area size to be searched；
Fig. 5 a is example diagram 1 corresponding when carrying out oriented search to Fig. 4；
Fig. 5 b is example diagram 2 corresponding when carrying out oriented search to Fig. 4；
Fig. 5 c is example diagram 3 corresponding when carrying out oriented search to Fig. 4；
Fig. 5 d is example diagram 4 corresponding when carrying out oriented search to Fig. 4.
Specific embodiment
To make the object, technical solutions and advantages of the present invention clearer, below in conjunction with attached in the embodiment of the present invention
Figure, is clearly and completely described technical solution of the present invention, it is clear that described embodiment is that a part of the invention is real
Example is applied, instead of all the embodiments.Based on the embodiments of the present invention, those of ordinary skill in the art are not making creation
Property labour under the premise of every other embodiment obtained, shall fall within the protection scope of the present invention.
As Figure 1-Figure 2, the Parallel districts implementation method of the present invention for SPH algorithm, includes the following steps,
(1), it creates the simulation model of pending simulation analysis and Initialize installation is carried out to the simulation model created；By
It can be compatible with finite element stimulation in the pre-treatment of SPH simulation calculation process, it can determine that it is initial by grid
Particle distribution situation, that is, need to complete the profile set of particle, the setting of simulated conditions, the attribute of material representated by particle and
The setting such as original state of particle.Specifically, as preferred embodiment of the invention, the step 1 comprising:
Step 101, creation simulation model simultaneously set grid file corresponding to created simulation model, obtain corresponding
Gridding information is simultaneously stored in interim grid data, and the gridding information includes grid serial number and grid material number, and grid
Information include with its corresponding to node number；
Node file corresponding to the created simulation model of step 102, setting, obtains corresponding nodal information and is stored in
In interim node data, the nodal information includes node number and node coordinate；
The centre coordinate of step 103, each grid in conjunction with corresponding to gridding information and nodal information phantom
And volume data corresponding to each grid, and corresponding equivalent diameter is calculated based on volume data calculated；
Above-mentioned each grid conversion is corresponding particle by step 104, wherein the position of each particle is i.e. corresponding to it
The centre coordinate of grid, the equivalent diameter of grid of the size i.e. corresponding to it, grid material of the material number i.e. corresponding to it
Number；
Step 105, setting material file, obtain the density of every kind of material corresponding to simulation model, Poisson's ratio, springform
The data such as the type and parameter of amount and material model；
Step 106 joins the attribute of set material and model according to the corresponding grid material number of each particle
Number is assigned to each particle；
Step 107, conditions setting and primary condition obtain condition value and conditioning corresponding to simulation model
Group piece number, there are in interim initial BVP data；
Step 108, assignment component message file simultaneously obtain package count corresponding to simulation model, if opening component in order
File then successively sets corresponding group of piece number, component type；Further, if current group piece number and component type satisfaction are worked as
Preceding corresponding first boundary condition then continues the member for setting current component, the serial number of related grid is obtained, by initial BVP condition
It is assigned to the corresponding particle of grid；If boundary condition acts on node, the cell data that the node belongs to is obtained, by boundary condition
It is assigned to the particle that corresponding unit is converted to；
Step 109, the initial control parameter of setting；The calculating control parameter include current time, it is overall calculate duration, when
Preceding time step, overall recurring number, current time step number, minimum time step-length, maximum time step-length, symmetric information and preservation knot
The information such as fruit setting, parallel information.
Step 1010, the attribute and model parameter for setting every kind of material, set material parameter and model parameter are pressed
Each particle is assigned to according to the material number of particle.
(2), the pressure value of each particle corresponding to the simulation model after calculating initialized setting, stress tensor and
Pseudo-viscosity value；Specifically, as preferred embodiment of the invention, the calculating process of the pressure data are as follows: first to each
Particle is recycled and is obtained the density and specific internal energy of each particle, is joined data obtained as the input for calculating pressure
Number；Secondly the type of its calculation of pressure function is determined in the state equation assigned according to each particle；Finally based on determined by
Parameter and calculation of pressure function are fully entered, terminates particle cyclic process after calculating the pressure of each particle；Specifically, conduct
Preferred embodiment of the invention, the calculating process of the stress value are as follows: each particle is recycled first and obtains each grain
Strain rate, specific rotation and the stress of son；Stress rate corresponding to each particle is calculated secondly based on stress rate formula；Again
After determining current time step, the stress of currently calculated particle is obtained multiplied by corresponding stress rate with time step
Amount, and by the way that the stress of the particle is obtained interim updated stress plus stress tensor, the stress of the particle is brought into
Set yield condition judges whether to surrender, the stress after being surrendered if surrender with radial return mapping method, if do not had
There is surrender then without processing；Finally the stress after judgement as the updated stress of the particle and is terminated particle and is circulated throughout
Journey；Specifically, as preferred embodiment of the invention, the calculating process of the pseudo-viscosity value are as follows: first to each effect into
Row recycles and successively selectes the particle i and particle j of one of effect centering, calculates the mutual work of this pair of effect particle
Firmly, that is, the relative coordinate (x of particle i and particle j is calculatedi-xj, yi-yj, zi-zj) and relative velocity (vxi-vxj, vyi-vyj,
vzi-vzj), wherein xi, yiAnd ziIndicate the direction x of particle i, the direction y and the direction z coordinate, xj, xjAnd zjIndicate the side x of particle j
To, the direction y and the direction z coordinate, vxi, vyiAnd vziIndicate the direction x of particle i, the direction y and the direction z speed, vxj, vyjAnd vzjTable
Show the direction x of particle j, the direction y and the direction z speed；Then calculate the smooth length h of particle ii, density rhoi, velocity of sound ciWith j
The smooth length h of sonj, density rhoj, velocity of sound cj；The mutual active force of this pair of effect particle is finally calculated, and should
Pseudo-viscosity power of the active force as this effect pair, and end effect is to circulation.
(3), in each set time step, the simulation model is scanned for analyzing, by determination it is each in terms of
Calculate neighbour's particle corresponding to particle；The searching analysis includes that step 31 determines each Parallel districts corresponding to simulation model
And each Parallel districts boundary；
Specifically, the step 31 includes the following steps: as preferred embodiment of the invention
Step 311 first recycles all particles corresponding to simulation model, finds the X in three-dimensional space respectively
Direction, Y-direction, in Z-direction, particle II corresponding to particle I corresponding to coordinate minimum and coordinate maximum；Secondly by particle
The coordinate of II and the coordinate of particle I subtract each other, and obtain the maximum span value on tri- directions X, Y, Z；Finally by by the emulation
The coordinate of the corresponding each particle of model successively subtract each other with the coordinate of particle I and divided by the coordinate direction obtained most
Large span value determines the opposite dimensionless coordinate (X of each particlej-Xi/X0,Yj-Yi/Y0,Zj-Zi/Z0), wherein Xj、Yj、Zj
Respectively indicate X-direction, Y-direction, the particle coordinate in Z-direction；Xi、Yi、ZiRespectively indicate X-direction, Y-direction, coordinate in Z-direction
The coordinate of the smallest particle I, X0Indicate the maximum span value in X-direction, YjIndicate the maximum span value in Y-direction, Z0Indicate Z
Maximum span value on direction, j >=1；
Step 312 calculates separately square of the dimensionless coordinate in X-direction, Y-direction, Z-direction corresponding to all particles
With, and the vector direction that the smallest dimensionless coordinate of quadratic sum is constituted is determined as principal direction；
Step 313, based on identified principal direction, according to the relative coordinate of all particles and principal direction from as low as big progress
Sequence: further, as preference of the invention, the step 313 is using the bucket sort method of floating number to all particles
It is ranked up, it is following to be illustrated by taking specific example as an example: to obtain the opposite dimensionless of all particles in a main direction first and sit
Scale value will obtain after opposite dimensionless coordinate value carries out coordinate format and be converted to corresponding character string, due to each conversion
The length having the same of character string afterwards, therefore all opposite dimensionless coordinates will not be negative value, and each character string is suitable
It is 0 that sequence, which is the 1st, and the 2nd is decimal point, and the 3rd to the 8th is 6 bit values after floating number decimal point, and the 9th is e, i.e. table
Show that index starts, the 10th be index symbol, be negative value, the 11 to 12nd be index ten and a position；Secondly, setting
After m bucket (this example is arranged 0 to No. 9, altogether 9 buckets), first the decimal place of each character string is ranked up: successively small to each
When numerical digit recycles, each particle is put into the bucket of corresponding serial number, and first put the particle first searched by set sequence
Enter i.e. if recycled to the 1st decimal, first decimal place 4 of current particle is put into No. four buckets, completes first decimal
After the circulation of position, the particle in each bucket is sequentially exported according to 0 to 9 sequence and (that is being exported since No. 0 bucket, is first put
Into the first output of No. 0 bucket) to form a new particle sequence, then start to carry out next bit decimal place with this new sequence
Sequence, the sequence until completing all decimal places；Then started pair with the corresponding output sequence of aforementioned completion decimal place sequence
Exponent bits are ranked up, but sort method at this time should be with decimal on the contrary, this is because the symbol of exponent bits is negative sign, usually
This number of bigger exponential representation is smaller, so corresponding relationship when being sequentially put into each particle barrel to each exponent bits circulation
Are as follows: if currently the value of sequence position is n first, and particle is put into 9-n bucket；It then will be in each bucket according to the sequence being put into
Particle output, the initiation sequence that sequence obtained sorts as next exponent bits at this time, until completing all exponent bits rows
Sequence, sequence obtained is exactly the particle sequence of all particles from small to large according to principal direction relative coordinate at this time；
Step 314 determines Parallel districts quantity to be delimited and determines primary number corresponding to each Parallel districts
And boundary particle serial number corresponding to each Parallel districts boundary position；Further, as the preferred step of the invention
314 include: step 3141, obtain determine simulation model corresponding to Parallel districts quantity；Such as work as according to current simulated environment
The configuration parameter of preceding simulation computer (especially cpu) determines Parallel districts quantity；Step 3142, according to corresponding to simulation model
All populations and the identified Parallel districts quantity integer obtained that is divided by add corresponding to each Parallel districts of a setting
Initial population；Step 3143 determines boundary particle serial number corresponding to each Parallel districts boundary position, i.e., according to first
Possessed by first Parallel districts of boundary particle serial number corresponding to a Parallel districts and the second Parallel districts boundary position
Population, and the n-th -2 Parallel districts of boundary particle serial number in (n-1)th parallel area and n-th of parallel area with (n-1)th simultaneously
The boundary particle serial number of row subregion recycles all Parallel districts plus the strategy of the population in (n-1)th parallel area,
With boundary particle serial number corresponding to each Parallel districts boundary position of determination, wherein n >=3, such as second Parallel districts and third
First Parallel districts of boundary particle serial number of a Parallel districts and the boundary particle serial number of second Parallel districts add second
Population possessed by a parallel area；
Step 315 determines each Parallel districts boundary position；Specifically, as preferred embodiment of the invention, the step
315 include: to be traversed according to the sequence that step 313 sequences, when traversing the boundary particle serial number of some Parallel districts,
The number of particle corresponding to the boundary particle serial number is extracted, and finds the absolute seat of the corresponding particle of the number in a main direction
Mark, this absolute coordinate is the boundary position of present parallel subregion Yu next Parallel districts；Complete time of all particle sequences
It goes through to obtain the coboundary of all Parallel districts and lower boundary, wherein coboundary is that current Parallel districts are divided parallel with previous
The boundary coordinate in area, lower boundary are the boundary coordinate of itself and latter Parallel districts；And first Parallel districts coboundary is main
Particle coordinate corresponding to the initial position of direction, the last one Parallel districts lower boundary are grains corresponding to principal direction final position
Subcoordinate；
Due to the relativeness between each Parallel districts be it is specific, orderly, by upper at the beginning of calculating
After the method for stating determines each Parallel districts, only need to be ranked up the particle inside each Parallel districts in subsequent calculating,
Then the number of particles for communicating to determine each Parallel districts by each Parallel districts changes, and changes adjustment boundary position according to quantity
It sets, i.e., content described in the described step 316；
Step 316 delimited Parallel districts and be ranked up to the particle of each Parallel districts, and logical by each Parallel districts
Letter is adjusted simultaneously to determine the number of particles situation of change in each Parallel districts based on identified number of particles situation of change
Row partition boundaries position.Specifically, the step 316 includes: step 3161, delimit parallel as preferred embodiment of the invention
Subregion, and particle possessed by each Parallel districts is ranked up from as low as big strategy according to principal direction, it is current to determine
Parallel districts corresponding to particle sequence；Step 3162 determines present parallel subregion after time step iterative cycles each time
Interior number of particles situation of change: i.e. since the first Parallel districts, judge whether particle corresponding in current region is more than
Primary number is, under the particle having more in particle sequence is given next Parallel districts (i.e. the second Parallel districts) and updated
The coboundary of one Parallel districts and the lower boundary of present parallel subregion；Step 3163 repeats step 3162, each to be sequentially completed
Number of particles situation of change corresponding to Parallel districts, and based on identified number of particles situation of change adjustment Parallel districts
Lower boundary position, such as receives the population that the second Parallel districts are counted after the particle that the first Parallel districts mark off, and judgement is current
Whether corresponding particle is more than primary number in region, is to give down the particle having more in identified particle sequence
One Parallel districts (i.e. third Parallel districts) simultaneously update the coboundary of third Parallel districts and the lower boundary of the second Parallel districts, directly
To the last one Parallel districts, corresponding number of particles is only counted at this time.
Specifically, the step 3 further includes step 32 to each Parallel districts carry out area as preferred embodiment of the invention
Domain search specifically includes step 321, by particle corresponding to each Parallel districts, is divided into the Parallel districts internal particle
And partition boundaries communicate particle；Division rule is to traverse to particle corresponding to current line subregion, and more current
Whether particle is less than itself set compacted support to the coboundary of the affiliated Parallel districts of the particle or the distance of lower boundary, if its
Distance to coboundary is less than itself compacted support, then the particle is confirmed to be the partition boundaries communicated with a upper Parallel districts
Communicate particle；If it is less than itself compacted support to the distance of lower boundary, which is confirmed to be carries out with next Parallel districts
The partition boundaries of communication communicate particle, if current particle is to the coboundary of the affiliated Parallel districts of the particle, the distance of lower boundary
It is not less than itself set compacted support, then is confirmed to be the Parallel districts internal particle；Step 322, respectively to per together
The internal particle of row subregion scans for, corresponding search strategy are as follows: step 3221 determines meter corresponding to present parallel subregion
Calculate domain size；Identified computational domain is divided into several regions according to a certain percentage by step 3222；Step 3223, will be each
Region is divided into several identical background grid units, determines the direction of search and calculates each area based on the identified direction of search
Number of meshes corresponding to domain further reduces aforementioned if the number of meshes is greater than set upper limit threshold
Ratio simultaneously repeats step 3222；Step 3224 determines region of search and searches for and closes on particle corresponding to current region；Described search
The size in domain is greater than current area size to be searched, so that described search domain can be closed on current region to be searched
Region partly overlaps；Step 3225, circulating repetition step 3224 are until the search of whole region is completed in search, and by search process
Obtained in each particle information be updated；Step 3226, by the side of the lower boundary of the Parallel districts after the completion of current search
Boundary's particle information passes to next Parallel districts.
It should be understood that this part is scanned in the internal particle to each Parallel districts or referred to as calculate node
When, be also what subregion (i.e. region described in step 3222) carried out in the search on this node, this Combined With Area Searching with
The search of aforementioned Parallel districts is different, at this time we assume that the particle in preceding calculate node is exactly all particles, passes through setting
Region of search or to close on particle corresponding to window search domain search current region, and search window region of search is only needed every time
In particle, after successively moving window region of search search again for the particle in other windows, the benefit done so is exactly to be not required to
The bounding box of all particles in calculate node is all once drawn upper background grid unit, only need to have grain in the window every time
Scanned in the case where son, and if window search domain is more much larger than background grid unit, it is extremely uneven in particle distribution
The window search domain of a large amount of not particles can once skip in the case where even, save a large amount of search time, and each time
Background grid unit only divides in window search domain, and the size in window search domain is more much smaller than whole bounding box, then may be used
To save many memory spaces, in this way while guaranteeing has the certain advantages of tree search, and do not need to face recursive procedure
The middle stack overflow occurred or the too deep problem of recurrence；It is experimentally confirmed that optimal region of search quantity should in each region of search
The quantity of background grid is suitable.It is scanned for simultaneously in region of search, in order to enable can be searched between each region of search, it is required that
The content of other region of search of covering part, i.e. determination described in step 3224 are all repeated when each search of each region of search of region of search
Region of search is simultaneously searched for when closing on particle corresponding to current region, and it is big that the size in described search domain is greater than current region to be searched
It is small, so that described search domain can partly overlap with the close region in current region to be searched, example diagram as shown in Figure 3,
Computational domain is divided into 4 pieces of impartial regions and is divided into each region several when scanning for any one piece of region
The identical background grid unit of layer and make the region of search slightly larger than the region i.e. its to there is with other regions a part Chong Die
The region of intersection, the region that preferably may make the overlapping to intersect is one layer of background grid unit, as shown in figure 4, it is marked
The region of horizontal line is the size in current region to be searched itself, and the region for being labeled with oblique line is the area for being additionally required search
Domain, to guarantee each interregional realization search；In order to search for the problem of domain search is repeated there is no search, we borrow oriented search
Method, such as example as shown in Figure 5, when so that search being scanned for along direction shown in Fig. 5 a, if current search
It is same background grid unit (by taking that background grid at center as an example), then needs to judge the corresponding node of the background grid
Number size, (such as be all when searched node number is bigger be a pair of of effect according to unified size criteria judgement
It is right) avoid repeat search；(it is greater than current partition on the left of region of search when central node (i.e. destination node) is located at left border
Part) when searching route it is as shown in Figure 5 b, and the grid where itself is not searched for, to prevent repeat search；Work as center
When node is located below boundary, searching route such as Fig. 5 c does not search for grid where itself equally；When central node is located at angle
It falls in boundary, searching route such as Fig. 5 d, only one search grid, does not need comparison section with the search between grid known to upper example
Point number ensures that not repeat search, and many times and computing resource can be saved by doing so, and searches for compared to traditional chained list,
Speed will not be reduced, and the demand of memory largely reduces.
(4), state equation and governing equation corresponding to phantom, specifically, as of the invention preferred
Example, when being calculated, if simulation model is symmetry model there are symmetrical border, that is, simulation model, to simulation model
It is acted on and information supplement is handled；Described acted on does not need information supplement processing to carry out more searches, only needs true
Recognize corresponding mirror, specific:
It is described to simulation model acted on to information supplement processing include:
Step 41 independently confirms each Parallel districts with the presence or absence of a symmetrical border, if current Parallel districts
There are symmetrical borders, then to all effects corresponding in current Parallel districts to recycling, and judge a certain effect pair
In with the presence or absence of at least one particle to symmetrical border distance be less than itself compacted support, be then confirm there may be with mirror image
The particle pair of effect；
Step 42, the mirror for further determining that the effect pair only need to consider one because being symmetrical relationship
The mirror image of a particle whether be another particle effect pair, that is, judge the mirror image particle of any one particle of the effect centering
Whether be another particle effect pair, the determination strategy are as follows: obtain the mirror image grain of any one particle of the effect centering
The mirror position of son, and judge acquired mirror position and the spacing of another particle and the difference of the sum of two particle compacted supports
Absolute value whether be less than identified judgment threshold, be then determining effect to being the particle with mirror image effect to confirming
The mirror image particle of any one particle of the effect centering and another particle act on pair each other, and if the effect centering there are it to arrive
The distance of partition boundaries is less than itself compacted support, then confirms the mirror image particle of the particle and the particle acts on pair each other；It is successively right
All effects act on to information simulation model circulation to obtain whole effects about symmetrical border to information
Carry out augment.
Further, as a preference of the present invention, due to described symmetrical not there is only 1/2 to simulation model, and exist
Multiple symmetrical border, such as when that is it is 1/4 symmetrical that symmetrical border, which is two, the mirror position pair of first time generation
A mirror position can be equally generated in second plane of symmetry, it is also possible to be the position of an effect particle for this mirror position
It sets, in addition there are 1/8 symmetrical situations can carry out effect pair for the symmetrical border of above-mentioned three kinds of situations according to following processes
Information supplement processing；
Step 41 independently confirms each Parallel districts with the presence or absence of symmetrical border, if current Parallel districts exist
Symmetrical border and symmetrical border number is at least 2, then respectively along the direction x, the direction y and the direction z to current parallel point
Corresponding all effects are to recycling in area, and during successively recycling to above-mentioned 3 directions, judgement is current
Effect centering be respectively less than itself compacted support with the presence or absence of the distance of at least one particle to wherein any one symmetrical border, be
Then there may be the particles pair with mirror image effect for confirmation；
Step 42, the mirror for further determining that the effect pair need to consider one because being symmetrical relationship
Mirror image of the particle on each plane of symmetry whether be another particle effect pair, that is, judge any one of the effect centering
The mirror image particle on each plane of symmetry of particle whether be another particle effect pair, the determination strategy are as follows: first
The mirror position of any one particle mirror image particle on a certain plane of symmetry of the effect centering is obtained, and judges acquired mirror
Whether the absolute value of the difference of the spacing and the sum of two particle compacted supports of image position and another particle is less than identified judgement
Threshold value is that then determining confirmation particle mirror image particle and another particle on the current plane of symmetry acts on pair each other, and if the work
With centering, there are its distances to partition boundaries to be less than itself compacted support, then confirms mirror image of the particle on the current plane of symmetry
Particle and the particle act on pair each other；Next obtains any one particle of effect centering mirror image particle on other planes of symmetry
Mirror position, and judge the spacing of acquired mirror position and another particle and the difference of the sum of two particle compacted supports
Whether absolute value is less than identified judgment threshold, be then determining confirmation particle on the plane of symmetry mirror image particle with it is another
Particle acts on pair each other, and if the effect centering there are its distances to partition boundaries to be less than itself compacted support, confirm the grain
Mirror image particle of the son on the plane of symmetry and the particle also effect pair each other；It is complete to obtain successively to all effects to circulation
The effect about symmetrical border in portion is acted on augmenting to information to information, to simulation model；The effect simultaneously
Include the label and respectively corresponding action direction of the particle of interaction to information, and acquisition effect to information it
Afterwards, pseudo-viscosity and interaction function are calculated according to symmetrical mark, is equally recorded in and acts among to information.
It is above-mentioned that the treatment process of multiple symmetrical border can be indicated by following step and formula:
(1) direction x, the direction y, the direction z are recycled respectively, is confirmed whether that there are symmetrical borders；I.e.
The direction z is recycled: for i=0 to sz
The direction y is recycled: for j=0 to sy
The direction x is recycled: for k=0 to sx
(2) judge whether a series of a series of mirror image particle A points of the A points of effect centering on each plane of symmetry are a little
B(Xb, Yb, Zb) effect pair, wherein a series of A point be coordinate set, be represented by (sign (k) * Xa, sign (j) * Ya,
sign(i)*Za)；Wherein, the direction x, the direction y, the direction z symmetrical mark be respectively sx, sy, sz, the 0 no symmetrical border of expression, 1
Indicate symmetrical border；Sign is symmetry direction symbol and sign (0)=1, sign (1)=- 1；* it indicates to be multiplied；It is to record A
Point and corresponding symmetry direction symbol sign are that a pair acts on to and calculates this pair of corresponding kernel function information.
(3) terminate the direction x circulation；Terminate the direction y circulation；Terminate the direction z circulation.
Specifically, the step 4 further includes physics control corresponding to phantom as preferred embodiment of the invention
Equation processed, the calculating of the physical control equation includes calculating for the calculating of solid relevant portion and fluid relevant portion, right
It answers, the solid relevant portion calculating includes: a1, continuity equation calculating, such as first to all effects to circulation, to obtain
The coordinate of particle, speed and symmetrical symbols are acted on, by the coordinate and speed that symmetrically meet composite calulation particle；Then it obtains
This to effect pair interaction function dwdx；This is obtained to the two of effect pair multiplied by action function with the relative velocity after synthesis
The rate of change of the density increment of a particle；After the increment that each particle obtains in accumulation loop calculating；End loop, so that each parallel
By stages is communicated；B1, momentum conservation equation calculate: to all effects to circulation, obtaining the stress of each effect particle, symmetrically
Symbol passes through symmetrical symbols composite calulation particle stress；Then the function dwdx of each effect Thermodynamic parameters is obtained；With synthesis
The acceleration increment of first particle, is closed based on this increment by symmetrical symbols in every a pair of of the particle of Stress calculation afterwards
At the acceleration increment of second particle；So that accumulation loop calculates increment corresponding to acquisition respectively to each particle alone；Knot
Shu Xunhuan, end loop, so that the acceleration of all particles is communicated and synthesized between each Parallel districts；C1, conservation of energy side
Journey calculates: to all effects to circulation, acquisition acts on the stress of particle, speed, and symmetrical symbols pass through symmetrical symbols composite calulation
With particle stress and speed；Obtain the function dwdx of each effect Thermodynamic parameters；With after synthesis stress and speed calculate it is each
To the specific internal energy change rate increment of particle in particle；So that accumulation loop calculates the increment obtained to each particle alone；End follows
Ring, so that the specific internal energy change rate of all particles is communicated and synthesized between each Parallel districts；It is corresponding, the fluid dependent part
Dividing calculating includes: a2, momentum conservation equation calculating: to all effects to circulation, the pressure that acquisition acts on particle is acted on this
Pair pseudo-viscosity, symmetrical symbols, when not using pseudo-viscosity, the pressure for passing through two particles synthesizes the pressure of Riemann Solution；
Obtain the function dwdx of each effect Thermodynamic parameters；The acceleration increment of first particle in every a pair of of particle is successively calculated, with
The acceleration increment of second particle is synthesized based on this increment by symmetrical symbols；So that each particle accumulation loop alone
Calculate the increment obtained；End loop, so that the acceleration of all particles is communicated and synthesized between each Parallel districts；B2, energy
It measures conservation equation to calculate: to all effects to circulation, obtaining the pressure of each effect particle and the pseudo-viscosity of this effect pair,
Speed, symmetrical symbols are synthesized when not using pseudo-viscosity by pressure by symmetrical symbols composite calulation particle rapidity
The pressure of Riemann Solution；Obtain the function dwdx of each effect Thermodynamic parameters；The ratio of each pair of particle is calculated with the speed after synthesis
Change of internal energy rate increment；So that accumulation loop calculates the increment obtained to each particle alone；End loop, so that each Parallel districts
Between communicated and synthesized the specific internal energy change rates of all particles.
Specifically, the step 4 further includes by carrying out time integral with renewal speed as preferred embodiment of the invention
With coordinate and other states and physical quantity: the update such as to speed and coordinate: to all real particles of all Parallel districts into
Row circulation；Obtain the acceleration, present speed and coordinate of each particle；So that each updated speed of particle adds equal to present speed
The product of upper acceleration and time step, updated coordinate are equal to changing coordinates plus updated speed and time step
Product；Such as to the update of energy variation: to all real particle circulations of all Parallel districts, obtaining the current specific internal energy of each particle
With specific internal energy change rate so that update after each particle specific internal energy be equal to current specific internal energy add specific internal energy change rate and time step
Long product；Such as to the update of strain: being recycled to all real particles of all Parallel districts；Obtain that each particle is current to answer
Change and strain rate；Each components of strain are recycled, so that the updated components of strain are equal to the current components of strain plus correspondence
Strain rate and time step product；Such as to the update of failure state: being followed to all real particles of all Parallel districts
Ring obtains the set formula of the updated stress and strain substitution of particle and judges particle failure state, marks if failure
For the particle that fails；Such as to the update of erosion state: being recycled to all real particles of all Parallel districts, obtain each particle
Time step, the information such as strain, substitutes into corresponding formula and judges erosion state, is labeled as corroding grain if reaching erosion standard
Son recycles again and is deleted to all particles the particle for corroding label, is finally reentered into the particle not corroded continuously
Storing data structure in as next step calculate particle data.
Specifically, the step 4 further includes carrying out boundary condition to calculate i.e. according to set by as preferred embodiment of the invention
Fixed boundary condition, handle the particle to suffer restraints speed and carry out coordinate update, specifically include the following steps: to it is all simultaneously
All real particle circulations of row subregion；Judge whether current particle is the particle for being marked as boundary marker, if particle is not
The particle of boundary marker then normally updates the speed and coordinate of the particle；If particle is the particle of boundary marker, further
Judge whether the displacement of the particle reaches the constraint of boundary condition, if reaching constraint, the displacement of the particle and coordinate are no longer
It updates, and speed is set to 0, if not reaching constraint, the speed and coordinate of the particle normally update.
Specifically, being illustrated by taking the relevant propagation of explosion route searching of detonation as an example: according to direct search method:
The point initiation of the initial setting up of the model is then determined first；Then the time of burst time and current time are obtained
Difference；It obtains the explosion velocity of explosive and determines initiation conditions, the determination of the initiation conditions refers to if particle is small to fire point position
In or equal to explosion velocity multiplied by the time difference, then reaches initiation conditions, otherwise do not reach initiation conditions；All particles are finally traversed,
If current particle reaches initiation conditions, which is set for detonation particle, is otherwise provided as the particle that do not detonate.According to
It connects searching method: not needing initial point initiation due to searching for indirectly, be the relative position by particle and the particle that detonated
Relationship determines, therefore comprising: successively judges whether to the model in addition to mirror image acts on all effects to other than to circulation
A pair of effect has detonated particle and the particle that do not detonate in, continues under judgement if being unsatisfactory for aforementioned condition
One effect pair；The burst time that one of detonation particle is obtained if meeting condition, calculate the time with current time
Difference, with the time difference multiplied by explosion velocity obtain detonation wave propagate distance, and further judge this distance with this to effect between
Distance, if interparticle distance from detonation wave propagation distance is greater than, carries out lower a pair of of effect to judgement, otherwise do not detonate particle
It will be initiated, and record interim detonation mark (interim detonation mark will not influence the judgement searched for below)；Finally complete all works
It is so far completed using interim detonation mark more new particle detonation mark so that each Parallel districts are communicated with to after circulation
The search process of detonation propagation of explosion.
(5), the simulation model is post-processed, that is, is based on set termination condition, exports corresponding simulation result.
Specifically, as preferred embodiment of the invention, if the step 5 includes: when the time currently calculated meeting output condition,
Emulation data corresponding to all Parallel districts are merged under set host process and are obtained the serial number of current output file
Label；All particles in aforementioned host process are recycled；File is written into the variable of the corresponding output of each particle and name variable；
After completing the data that display needs, same text is written into the increment information for restarting needs and material and calculating control information
Part；Corresponding simulation result is finally exported, and updates the output serial number label of current file.Specifically, as of the invention excellent
Example is selected, the step 5 includes that setting print screen calculates progress, it may be assumed that is communicated to all Parallel districts, acquisition is each simultaneously
Print screen corresponding to row subregion needs the value of parameter；Each value obtained is screened, obtains maximum value therein, most
Small value and the average value being respectively worth；Each value after screening is substituted into the function of set print screen；According to set defeated
Line number out adjusts cursor position, the content for enabling refreshing just covering each time need to refresh, to prevent asking for wrong row
Topic；Execute print screen task.
The foregoing is only a preferred embodiment of the present invention, but scope of protection of the present invention is not limited thereto,
Anyone skilled in the art in the technical scope disclosed by the present invention, according to the technique and scheme of the present invention and its
Inventive concept is subject to equivalent substitution or change, should be covered by the protection scope of the present invention.
Claims (5)
1. a kind of Parallel districts implementation method for SPH algorithm, it is characterised in that include the following steps:
Step 1, the simulation model of the pending simulation analysis of creation simultaneously carry out Initialize installation to the simulation model created；
Step 101, creation simulation model simultaneously set grid file corresponding to created simulation model, obtain corresponding grid
Information is simultaneously stored in interim grid data, and the gridding information includes grid serial number and grid material number, and gridding information
Comprising with its corresponding to node number；
Node file corresponding to the created simulation model of step 102, setting obtains corresponding nodal information and is stored in interim
Node data in, the nodal information includes node number and node coordinate；
Step 103, each grid in conjunction with corresponding to gridding information and nodal information phantom centre coordinate and
Volume data corresponding to each grid, and corresponding equivalent diameter is calculated based on volume data calculated；
Above-mentioned each grid conversion is corresponding particle by step 104, wherein grid of the position of each particle i.e. corresponding to it
Centre coordinate, the equivalent diameter of grid of the size i.e. corresponding to it, grid material number of the material number i.e. corresponding to it；
Step 105, setting material file, obtain simulation model corresponding to every kind of material density, Poisson's ratio, elasticity modulus with
And the type and supplemental characteristic of material model；
Step 106 assigns the attribute of set material and model parameter according to the corresponding grid material number of each particle
To each particle；
Step 107, conditions setting and primary condition obtain the component of condition value and conditioning corresponding to simulation model
Number, there are in interim initial BVP data；
Step 108, assignment component message file simultaneously obtain package count corresponding to simulation model, if opening component text in order
Part then successively sets corresponding group of piece number, component type；Further, if current group piece number and component type meet currently
Corresponding first boundary condition then continues the member for setting current component, obtains the serial number of related grid, initial BVP condition is assigned
Give grid corresponding particle；If boundary condition acts on node, the cell data that the node belongs to is obtained, boundary condition is assigned
The particle being converted to corresponding unit；
Step 109, the initial control parameter of setting；When the initial control parameter includes current time, totally calculates duration, is current
Between step-length, overall recurring number, current time step number, minimum time step-length, maximum time step-length, symmetric information and save result and set
Fixed, parallel information；
Step 1010, the attribute and model parameter for setting every kind of material, by set material parameter and model parameter according to grain
The material number of son is assigned to each particle；
Pressure value, stress tensor and the people of each particle corresponding to simulation model after step 2, the initialized setting of calculating
Work viscosity value；
Step 3 scans for analyzing in each set time step to the simulation model；The searching analysis packet
Include each Parallel districts corresponding to determining simulation model and each Parallel districts boundary；
Step 31 determines each Parallel districts and each Parallel districts boundary corresponding to simulation model comprising:
Step 311 first recycles all particles corresponding to simulation model, find respectively three-dimensional space i.e. X-direction,
In Y-direction, Z-direction, coordinate smallest particles I and coordinate maximum particle II；Secondly by the coordinate of above-mentioned particle II with particle I's
Coordinate subtracts each other, and obtains the maximum span value on tri- directions X, Y, Z；Finally by by each grain corresponding to the simulation model
The coordinate of son successively subtracts each other with the coordinate of particle I and divided by the maximum span value on coordinate direction obtained, determines and stores
The opposite dimensionless coordinate (X of each particlej-Xi/X0,Yj-Yi/Y0,Zj-Zi/Z0), wherein Xj、Yj、ZjRespectively indicate X-direction, Y
Particle coordinate on direction, Z-direction；Xi、Yi、ZiRespectively indicate X-direction, Y-direction, in Z-direction the smallest particle I of coordinate seat
Mark, X0Indicate the maximum span value in X-direction, Y0Indicate the maximum span value in Y-direction, Z0Indicate the maximum span in Z-direction
Value, j >=1；
Step 312 calculates separately quadratic sum of the dimensionless coordinate in X-direction, Y-direction, Z-direction corresponding to all particles,
And the vector direction that the smallest dimensionless coordinate of quadratic sum is constituted is determined as principal direction；
Step 313, based on identified principal direction, the relative coordinate of each particle in a main direction is found out, then according to opposite
Coordinate is ranked up all particles from as low as big sequence；
Step 314 determines primary number corresponding to each Parallel districts and each based on Parallel districts quantity to be delimited
Boundary particle serial number corresponding to the boundary position of Parallel districts；
Step 315 determines each Parallel districts boundary position and stores the corresponding boundary coordinate of identified boundary particle serial number；
Step 316 delimited Parallel districts and be ranked up to the particle inside each Parallel districts, and logical by each Parallel districts
Letter is adjusted simultaneously to determine the number of particles situation of change in each Parallel districts based on identified number of particles situation of change
Row partition boundaries position, to reach load balancing；
State equation corresponding to step 4, phantom and governing equation；
Step 5 exports simulation result based on set termination condition.
2. according to the method described in claim 1, it is characterized by:
All particles are ranked up using the bucket sort method of floating number in the step 313.
3. according to the method described in claim 1, it is characterized by:
The step 314 includes:
Step 3141 obtains Parallel districts quantity corresponding to simulation model；
Step 3142, be divided by according to all populations corresponding to simulation model with identified Parallel districts quantity it is obtained
Integer adds one to determine population initial inside each Parallel districts；
Step 3143 determines boundary particle serial number corresponding to each Parallel districts boundary position, i.e., according to first Parallel districts
With population possessed by first Parallel districts of boundary particle serial number corresponding to the second Parallel districts boundary position, and
Point in n-1 parallel areas and boundary the n-th -2 Parallel districts of particle serial number and (n-1)th Parallel districts in n-th of parallel area
Boundary's particle serial number recycles all Parallel districts plus the strategy of the population in (n-1)th parallel area, to determine respectively simultaneously
Boundary particle serial number corresponding to row partition boundaries position, wherein n >=3.
4. according to the method described in claim 1, it is characterized by:
The step 315 includes:
It is traversed according to the sequence that step 313 sequences, when traversing the boundary particle serial number of some Parallel districts, is extracted
The number of particle corresponding to the boundary particle serial number, and the relative coordinate of the corresponding particle of the number in a main direction is found,
This relative coordinate is the boundary position of present parallel subregion Yu next Parallel districts；Complete the traversal of all particle sequences with
Coboundary and the lower boundary of all Parallel districts are obtained, wherein coboundary is current Parallel districts and previous Parallel districts
Boundary coordinate, lower boundary are the boundary coordinate of itself and latter Parallel districts；And first Parallel districts coboundary is principal direction
Particle coordinate corresponding to initial position, the last one Parallel districts lower boundary are that particle corresponding to principal direction final position is sat
Mark.
5. according to the method described in claim 1, it is characterized by:
The step 316 includes:
Step 3161 delimit Parallel districts, and according to the certainly as low as big strategy of principal direction to grain possessed by each Parallel districts
Son is ranked up, with particle sequence corresponding to the current Parallel districts of determination；
Step 3162, the number of particles situation of change after time step iterative cycles each time in determining present parallel subregion: i.e.
Since the first Parallel districts, judge whether particle corresponding in current region is more than primary number, is then from particle sequence
The particle having more is given into next Parallel districts in column, and is updated under the coboundary and present parallel subregion of next Parallel districts
Boundary；
Step 3163 repeats step 3162, to be sequentially completed number of particles situation of change corresponding to each Parallel districts, and base
Parallel districts up-and-down boundary position is adjusted in identified number of particles situation of change.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610968969.7A CN106529011B (en) | 2016-11-03 | 2016-11-03 | A kind of Parallel districts implementation method for SPH algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610968969.7A CN106529011B (en) | 2016-11-03 | 2016-11-03 | A kind of Parallel districts implementation method for SPH algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106529011A CN106529011A (en) | 2017-03-22 |
| CN106529011B true CN106529011B (en) | 2019-05-10 |
Family
ID=58325891
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610968969.7A Active CN106529011B (en) | 2016-11-03 | 2016-11-03 | A kind of Parallel districts implementation method for SPH algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106529011B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109711525A (en) * | 2018-12-12 | 2019-05-03 | 湖北航天技术研究院总体设计所 | A kind of proximate particle search method and system for SPH algorithm |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7948485B1 (en) * | 2005-12-12 | 2011-05-24 | Sony Computer Entertainment Inc. | Real-time computer simulation of water surfaces |
| CN102262689A (en) * | 2010-05-26 | 2011-11-30 | 利弗莫尔软件技术公司 | Hybrid element enabling finite element/smoothed particle hydrodynamics coupling |
| CN102831280A (en) * | 2012-09-10 | 2012-12-19 | 北京航空航天大学 | Meshless physical deformation simulation method based on moving least squares |
| CN102930087A (en) * | 2012-10-19 | 2013-02-13 | 湖南大学 | Method for searching adjacent particles in analog simulation technology |
| CN103559741A (en) * | 2013-11-25 | 2014-02-05 | 武汉大学 | Particle-based multiphase coupling method in virtual surgery |
| CN104951617A (en) * | 2015-06-01 | 2015-09-30 | 电子科技大学 | Antenna optimization design method based on area decomposition method |
-
2016
- 2016-11-03 CN CN201610968969.7A patent/CN106529011B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7948485B1 (en) * | 2005-12-12 | 2011-05-24 | Sony Computer Entertainment Inc. | Real-time computer simulation of water surfaces |
| CN102262689A (en) * | 2010-05-26 | 2011-11-30 | 利弗莫尔软件技术公司 | Hybrid element enabling finite element/smoothed particle hydrodynamics coupling |
| CN102831280A (en) * | 2012-09-10 | 2012-12-19 | 北京航空航天大学 | Meshless physical deformation simulation method based on moving least squares |
| CN102930087A (en) * | 2012-10-19 | 2013-02-13 | 湖南大学 | Method for searching adjacent particles in analog simulation technology |
| CN103559741A (en) * | 2013-11-25 | 2014-02-05 | 武汉大学 | Particle-based multiphase coupling method in virtual surgery |
| CN104951617A (en) * | 2015-06-01 | 2015-09-30 | 电子科技大学 | Antenna optimization design method based on area decomposition method |
Non-Patent Citations (3)
| Title |
|---|
| An improved parallel SPH approach to solve 3D transient generalized Newtonian free surface flows;Jinlian Ren;《Computer Physics Communications》;20160428;第87-105页 |
| Towards consistence and convergence of conservative SPH approximations;S.Litvinov;《Journal of Computational Physics》;20151115;第394-401页 |
| 高速碰撞数值计算中的SPH分区算法;卞梁;《计算物理》;20110331;第28卷(第2期);第207-212页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106529011A (en) | 2017-03-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106528989B (en) | A kind of distributed parallel SPH emulation mode | |
| CN106503365B (en) | A kind of sector search method for SPH algorithm | |
| US11675940B2 (en) | Generating integrated circuit floorplans using neural networks | |
| CN106485030B (en) | A kind of symmetrical border processing method for SPH algorithm | |
| CN108122032A (en) | A kind of neural network model training method, device, chip and system | |
| JP2023522567A (en) | Generation of integrated circuit layouts using neural networks | |
| CN110135584A (en) | Extensive Symbolic Regression method and system based on self-adaptive parallel genetic algorithm | |
| CN106250933A (en) | Method, system and the FPGA processor of data clusters based on FPGA | |
| CN115437795B (en) | Video memory recalculation optimization method and system for heterogeneous GPU cluster load perception | |
| CN102521854A (en) | Parallel flow line placing method applicable to two-dimensional flow field | |
| CN112541584B (en) | Deep neural network model parallel mode selection method | |
| CN107229234A (en) | The distributed libray system and method for Aviation electronic data | |
| CN113283186A (en) | Universal grid self-adaption method for CFD | |
| CN109598742A (en) | A kind of method for tracking target and system based on SSD algorithm | |
| KR20220134627A (en) | Hardware-optimized neural architecture discovery | |
| CN111767689A (en) | Three-dimensional integrated circuit layout method based on graphic processing | |
| CN113435128A (en) | Oil and gas reservoir yield prediction method and device based on condition generation type countermeasure network | |
| Liu et al. | Quantum-inspired African vultures optimization algorithm with elite mutation strategy for production scheduling problems | |
| CN106529011B (en) | A kind of Parallel districts implementation method for SPH algorithm | |
| CN109961129A (en) | A kind of Ocean stationary targets search scheme generation method based on improvement population | |
| CN104698838B (en) | Based on the fuzzy scheduling rule digging method that domain dynamic is divided and learnt | |
| CN106484532A (en) | GPGPU parallel calculating method towards SPH fluid simulation | |
| CN109032667A (en) | Adjacency list method for fast establishing and system in a kind of molecular dynamics simulation | |
| Kampolis et al. | Multilevel optimization strategies based on metamodel-assisted evolutionary algorithms, for computationally expensive problems | |
| CN116339973A (en) | Digital twin cloud platform computing resource scheduling method based on particle swarm optimization algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |