JP2012128793A - Particle state calculation device and particle state calculation method - Google Patents

Abstract

PROBLEM TO BE SOLVED: To efficiently calculate a particle state through decentralized memory type parallel calculation.SOLUTION: A plurality of nodes 10 are apparatuses that constitute a particle state calculation system 1, and calculates states of particles arranged in their own regions as a plurality of divided calculation regions, in each time step. Each of the nodes 10 includes a particle state calculation unit 11 which calculates a state of particles in its own region; a period determination unit 12 which determines whether a time step is a time step having a predetermined period; an own-virtual-particle information output unit 13 which outputs information representing own region virtual particles used by other nodes in a time step of the predetermined period; an other-virtual-particle information input unit 14 which inputs other-region virtual particle information used for calculation of the own region; an own-virtual-particle state output unit 15 which outputs a state of own-region virtual particles; and an other-virtual-particle state input unit 16 which inputs states of other-region virtual particles. The particle state calculation unit 11 calculates the state of particles using information on other-area virtual particles.

Description

  The present invention relates to a particle state calculation apparatus and a particle state calculation method for calculating a particle state by distributed memory type parallel calculation.

  In the early days of cloud development, the particle system distribution of cloud particles (water particles) spreads rapidly (dispersion of the particle system distribution increases rapidly), and a phenomenon that causes a quick start of rainfall has been discussed for a long time. The causes include the effect of cloud particle collision acceleration by cloud turbulence, turbulent entrainment, giant condensation nuclei, and turbulent diffusion of cloud particles. Above all, the collision promotion effect by cloud turbulence mentioned at the beginning has been discussed for the longest time. Along with this, research on the collision frequency of particles in turbulent flow has been actively conducted.

  Some models have been proposed to predict the turbulent collision frequency of particles. For example, a method based on collision frequency factor data obtained by direct numerical simulation (DNS) as described in Non-Patent Document 1 has been proposed.

Woittiez, E.J.P., Jonker, H.J.J. and Portela, L.M., “On the Combined Effects of Turbulence and Gravity on Droplet Collisions in Clouds: a Numerical Study”, Journal of the AtmosphericSciences, Vol. 66 (2009), pp. 1926-1943.

However, the prediction results obtained by the above conventional method have not been confirmed to be reliable at high Reynolds numbers as follows. As shown in Non-Patent Document 1, the data obtained in the conventional research is a Taylor microscale reference Reynolds number Re _λ (= u′l _λ / ν, which is at most 85, where u ′ is the velocity fluctuation intensity, l _(λ is the Taylor microscale, ν is the kinematic viscosity coefficient). The Reynolds number is far from the high Reynolds number of Re _λ = 10 ^3-4 in actual cloud turbulence. Therefore, the validity of applying the conventional turbulent collision model to the cloud turbulence simulation as it is remains doubtful.

  In order to solve this, collision frequency factor data at a Reynolds number as high as possible is required. For this purpose, large-scale parallel computation is essential. Parallel computation includes shared memory parallel computation (openMP, automatic parallelization library, etc.) where each process shares the same physical memory, and distributed memory parallel computation (Message Passing Interface) where each process uses separate physical memory. (MPI) library is generally).

  The former can only execute parallel computation using tens of processes at most, and the memory that can be used is limited. Therefore, there is a disadvantage that large-scale parallel computation cannot be executed. So far, shared memory (openMP, etc.) parallel code for particle collisions in turbulent flow has been developed. However, for the above reasons, it is difficult to obtain particle collision data at a high Reynolds number in shared memory parallel computation.

  On the other hand, in the latter distributed memory type parallel computation, if hardware is gathered, parallel computation using many processes is possible. In the distributed memory type parallel calculation, each process uses a (distributed) memory, and advances the calculation while performing data communication when necessary. Since the data communication speed is extremely slower than the memory access speed, efficient data communication is required.

  However, a distributed memory type parallel computation code for particle collision has not been developed. One of the reasons is that DNS using a spectral method has been performed mainly for gas phase turbulent flow calculation. When spectrum calculation is performed in a distributed memory type parallel environment, communication of all area data (all-to-all) occurs. Since a technique for efficiently processing a huge amount of data communication is required, a spectrum method corresponding to three-dimensional region division has not been developed at present. In a large-scale parallel calculation based on two-dimensional region division, each region has a long, so-called pencil shape, and (surface area) / (volume) increases. In the case of parallel computation that requires passing information between neighboring processes such as particle computation, the amount of communication is proportional to the surface area of the region, so the increase in (surface area) / (volume) is (communication amount) / (computation amount). ), Resulting in a decrease in parallel efficiency. Therefore, two-dimensional division is not suitable for large-scale parallel computation for particle motion. From these facts, it is considered that there is a great difficulty in calculating turbulent flow using a spectral method that supports only two-dimensional division and acquiring particle collision data at a high Reynolds number.

  The present invention has been made in view of the above, and a particle state calculation apparatus and a particle state calculation that can efficiently calculate a particle state by distributed memory type parallel calculation and can perform calculation for a larger model than the result. It aims to provide a method.

  In order to achieve the above object, the particle state calculation apparatus according to the present invention is a particle state calculation that repeatedly calculates the state of particles arranged in its own region, which is one of a plurality of divided calculation regions, at each time step. A device for calculating a state of a particle in its own region at each time step, and a period for determining whether or not the target time step in the particle state calculation unit is a time step of a preset period Based on the position of the particle calculated by the particle state calculation means, the time step determined by the determination means and the time step of the period preset by the period determination means is adjacent to the own region from the particle. The self-region virtual particles used for the calculation in the adjacent region that is the divided calculation region are identified, and information indicating the identified self-region virtual particles is stored in the adjacent region. Self-virtual particle information output means for outputting to another particle state calculation device for calculating the state of the other, and other virtual for inputting information specifying other region virtual particles used for calculation in the own region from another particle state calculation device Particle information input means, self-virtual particle state output means for outputting information indicating the state of the self-region virtual particles calculated by the particle state calculation means at each time step to another particle state calculation device, and separate at each time step Particle state calculation means, comprising: other virtual particle state input means for inputting information indicating the state of the virtual particle in the other region from the particle state calculation apparatus, and output means for outputting information according to the calculation by the particle state calculation means Are information specifying other region virtual particles input by other virtual particle information input means, and other region virtual particles input by other virtual particle state input means And calculating the state of particles using the information indicating the state.

  In the particle state calculation apparatus according to the present invention, the local region virtual particles in the local region are specified in a time step of a preset period, and the information indicating the specified local region virtual particles calculates the state of the adjacent region. Is output to the particle state calculation device. On the other hand, information specifying other region virtual particles used for calculation in the own region is input from another particle state calculation device. In each time step, information indicating the state of the own region virtual particles is output and information indicating the state of the other region virtual particles is input. The state of the particle in its own region is calculated using the information specifying the other region virtual particles and the information indicating the state of the particle. That is, by making particles related to information exchanged between particle state calculation devices limited virtual particles, it is possible to reliably input necessary information in each particle state calculation device and efficiently communicate between particle state calculation devices. To. Further, communication between the particle state calculation devices becomes more efficient by exchanging information for specifying virtual particles as a time step having a preset period. Thereby, according to the particle state calculation apparatus according to the present invention, the particle state can be efficiently calculated by the distributed memory type parallel calculation, and the result can be calculated for a larger model.

  The particles are particles in the fluid, and the particle state calculation means preferably calculates the state of the fluid by a finite difference method and calculates the state of the particles in consideration of the influence of the fluid. According to this configuration, a large-scale model calculation related to particles in a fluid can be performed.

  The particle state calculation means specifies the proximity relation between particles based on the information specifying the other region virtual particles input by the other virtual particle information input means, and calculates the state of the particles based on the neighborhood relation. It is desirable. According to this configuration, it is possible to reliably perform calculation of a model in which interaction between particles such as collision between particles occurs.

  By the way, the present invention can be described as an invention of a particle state calculation apparatus as described above, and can also be described as an invention of a particle state calculation method as follows. This is substantially the same invention only in different categories, and has the same operations and effects.

  That is, the particle state calculation method according to the present invention is a particle state calculation method by a particle state calculation device that repeatedly calculates the state of particles arranged in its own region, which is one of a plurality of divided calculation regions, for each time step. A particle state calculation step for calculating the state of the particle in its own region at each time step, and a period determination for determining whether or not the target time step in the particle state calculation step is a time step of a preset period. A step adjacent to its own region from the particle based on the position of the particle calculated in the particle state calculation step to the time step determined to be a time step of a period set in advance in the step and the cycle determination step The self-region virtual particles used for the calculation are specified in the adjacent region that is the calculated region, and the identified self-region virtual particles are indicated. The self-virtual particle information output step that outputs the information to another particle state calculation device that calculates the state of particles in the adjacent region, and the other region virtual particles used for the calculation in the own region are identified from another particle state calculation device Another virtual particle information input step for inputting information to be performed, and a self virtual particle state output step for outputting information indicating the state of the own region virtual particles calculated in the particle state calculation step in each time step to another particle state calculation device And another virtual particle state input step for inputting information indicating the state of the other region virtual particles from another particle state calculation device at each time step, and an output step for outputting information according to the calculation in the particle state calculation step, In the particle state calculation step, other region virtual particles input in the other virtual particle information input step are specified. Information, and characterized by computing the condition of the particles using the information indicating the state of the input other areas virtual particles in the other virtual particle state input step.

  According to the present invention, the particles related to information exchanged between the particle state calculation devices are defined as limited virtual particles, so that necessary information can be surely input in each particle state calculation device and between the particle state calculation devices. Make communication efficient. Further, communication between the particle state calculation devices becomes more efficient by exchanging information for specifying virtual particles as a time step having a preset period. Thus, according to the present invention, it is possible to efficiently calculate a particle state by distributed memory type parallel calculation, and to calculate a model larger than the result.

It is a figure which shows the function structure of the node which is the particle | grain state calculation apparatus which concerns on embodiment of this invention. It is a figure which shows the division | segmentation of the calculation area | region (in fluid calculation) of a fluid and a particle | grain. It is a figure which shows the division | segmentation into the cell of a calculation area. It is a flowchart which shows the process performed with the whole particle state calculation system comprised by a some node. It is a flowchart which shows the process (particle state calculation method) performed by each node which is the particle state calculation apparatus which concerns on embodiment of this invention. It is a figure which shows the particle | grains in the calculation area | region according to each process of the flowchart of FIG. 5, and is made into the top view for description. It is a graph which shows the result of the collision frequency factor of the particle | grains obtained by the particle | grain state calculation apparatus and particle state calculation method which concern on this embodiment. It is a graph which shows the result of the collision frequency factor of the particle | grains whose Stokes number is 1 obtained by the particle | grain state calculation apparatus and particle | grain state calculation method which concern on this embodiment. It is a graph which shows the calculation time required for (a) 1 hour step at the time of changing the number of calculation processes with respect to N256 calculation, and (b) calculation speed.

  Hereinafter, preferred embodiments of a particle state calculation apparatus and a particle state calculation method according to the present invention will be described in detail with reference to the drawings. In the description of the drawings, the same elements are denoted by the same reference numerals, and redundant description is omitted.

  FIG. 1 shows a node 10 that is a particle state calculation apparatus according to the present embodiment. The node 10 constitutes the particle state calculation system 1 together with another node 10. The particle state calculation system 1 is a system that calculates (simulates) the state of particles such as water particles constituting a cloud, for example. Each node 10 is a computer having hardware such as a CPU (Central Processing Unit), a memory, and a hard disk. Each node 10 is connected by a network such as an intranet or the Internet, and information can be transmitted and received between the nodes 10. The calculation of the particle state is performed by a distributed memory type parallel calculation in which each node 10 performs a calculation by using a separate physical memory. Here, the calculation process in each node 10 is referred to as an (arithmetic) process.

  The particles to be calculated are arranged in a calculation area whose size is set in advance. The calculation area is divided into a plurality of areas. Each node 10 is assigned one of the divided areas as a calculation target, and repeatedly calculates the state of particles arranged in its own area, which is the assigned area, for each time step. The calculated particle state is specifically the position and velocity of the particle. In addition, a fluid (more specifically, an incompressible fluid) flows in the calculation region, and the particles may be in the fluid and affected by the fluid (that is, the fluid is also calculated). The fluid is, for example, an air flow that affects the state of water particles. The calculated fluid state is specifically the velocity and pressure of the fluid.

  The nodes 10 adjacent to each other in the calculation target area transmit and receive particle information according to the movement of the particle and the like, and use the particle information for the particle calculation.

  The particle state calculation system 1 also includes an apparatus for comprehensively controlling the distributed memory type parallel calculation. The apparatus performs processing for overall control of distributed memory type parallel calculation such as input of parameters for calculation, division of a calculation area, and calculation instruction to each node 10. The device may be any one of the plurality of nodes 10.

Here, calculation of particles and fluid will be described. First, the governing equation of the fluid and the numerical calculation method will be described. The incompressible fluid equation is the following continuity equation and the Navier-Stokes equation:

Here, u _i is a fluid velocity vector, and p is a variable indicating the fluid pressure. t is a variable indicating a time step. The subscripts i and j represent any direction of x, y, and z, and the expressions (1) and (2) are expressed according to the sum rules. All variables are made dimensionless by the representative length L ₀ and the representative speed U ₀ . Re is a value given as a parameter defined as Re = U ₀ L ₀ / ν using the kinematic viscosity coefficient ν. For example, ν is set to a value of 1.5 × 10 ⁻⁵ m ² / s at 1 atmosphere 298K. The last term f _{i on} the right side of Equation (2) is an external force term.

For the calculation of the convection difference, Morinishi, Y., Lund, T., Vasilyev, O. and Moin, P., “Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow”, Journal of Computational Physics, Vol. 143 ( 1998), pp. 90-124. (4th-order scheme for Non-linear term), a conservative fourth-order central difference scheme, and the fourth-order central difference method can be used for the calculation of the difference of the viscosity term. For the calculation of time progression (time evolution), the second-order Runge-Kutta method is used, and the HSMAC method (Hirt) is used for iterative calculation while simultaneously relaxing the velocity and pressure so that the velocity after time evolution satisfies the continuous equation. , CW and Cook, JL, “Calculating Three-Dimensional Flowsaround Structures and over Rough Terrain”, Journal of Computational Physics, Vol. 10 (1972), pp.324-340). As the spatial difference method in the iterative calculation, a second-order center difference method can be used. The convergence determination condition for iterative calculation can be, for example, ∂ _i u _i / (U ₀ / Δ) <10 ⁻³ (Δ is the calculation grid width).

The external force term f _i, can be used Reduced-CommunicationForcing (RCF) method (Onishi, R, Baba, Y. andTakahashi, K., "Efficient Large-Scale Forcing in Finite-Difference Simulationsof Steady Isotropic Turbulence", Topics on Chaotic Systems -selected papersfrom CHAOS 2010 international conference, Skiadas, CH et al. eds., World Scientific Publishing Co. (in press)). This is one of the large-scale forcing methods in which energy is injected only into large-scale motion, and is a forcing method that can achieve high parallel efficiency in the finite difference method. For example, the magnitude of the wave number k is set so that the energy of large-scale motion with | k | <2.5 is maintained.

Next, the fluid calculation area and the parallel calculation method will be described. As shown in FIG. 2, N ³ staggered lattices are arranged in a cubic calculation region having a side of 2πL ₀ . That is, the calculation lattice width Δ = 2πL ₀ / N. Here, N is a parameter set or input in advance. A specific calculation code can be described by a Fortran 90 code using an MPI library. This corresponds to three-dimensional area division, and the x, y, and z directions, which are coordinate axes of the cube calculation area, can be divided at equal intervals. The division number _M x in each direction, _{M y} and _{M z} are in relationships of MPI processes total _{M proc} and _{M proc = M x M y M} z. In this case, each process having an _{n x} × _{n y} × _n magnitude of _{z (n x = N / M} x, n y = N / M y and _{n z = N / M z)} "MPI divided region" Responsible for the calculation of fluid and particle motion within M _x , M _y, and M _z are parameters that are set or input in advance as arbitrary divisor values of N. In the above example, N ³ is divided at equal intervals. However, the number of divisions can be changed in each direction such as N _x , N _y , and N _z .

For example, the gas phase turbulent flow field is calculated under the calculation conditions shown in the following table. In all cases, the time increment Δt was set so that the Courant number (U ₀ Δt / Δ) was 0.2 or less. Here, the time step Δt is the time step interval of the repeated calculation. However, depending on the size of the particles, it may be necessary to set Δt to be smaller as dictated.

Subsequently, the calculation of the particle state will be described. Particle motion can be calculated, for example, using Lagrangian particle tracking. Assuming linear Stokes drag, the governing equation is

Here, v _p is a particle velocity, and u _i (the description of the underbar is omitted) is a variable indicating the fluid velocity at the particle position. Here, the subscript i represents the x, y, and z directions. τ _p is the particle relaxation time. Similar to equation (1), all variables are made dimensionless by the representative length L ₀ and the representative speed U ₀ . Also, the gravity acting on the particles is ignored. The particle relaxation time τ _p is calculated by the following equation.

Here, r is the particle radius, and ρ _p and ρ _f are the density of the particle and the fluid, respectively, which are parameters set or input in advance. For example, assuming a cloud particle, ρ _p / ρ _f = 840 from the density values of water and air at 1 atmosphere 298K. The linear Stokes drag model is applicable when the particle Reynolds number Re _p (= 2r | v _p −u | / ν) is sufficiently smaller than 1. Assuming cloud droplets, Re _p = 1 when r = 40 μm. For larger particles, a nonlinear drag model should be used, but for simplicity, a linear Stokes drag model may be used in all cases. Moreover, since the maximum mass ratio of water droplets to air is about 10 ^{−3 in} an actual cloud, the particles may be dilute and do not affect the fluid motion. The fluid velocity u at the particle position is calculated by linear interpolation from the eight velocity data in the vicinity of the particle. For the calculation of time progression, the second-order Runge-Kutta method is used as in the fluid calculation.

Next, particle collision calculation will be described. The collision considers only the collision between two particles. A trajectory in which two particles move between Δt is obtained, and it is determined whether or not the minimum distance of these trajectories is equal to or smaller than the particle diameter. One particle disappears after the collision. Therefore, once a collision occurs, the total number of particles decreases by one. Using this, the collision frequency N _c and the collision frequency factor K _c (= 2N _c / n _p ² , n _p is the particle number density) are calculated from the decrease rate of the total number of particles. From the particle total at a certain time step n _N ^p n, 1 Step particles total number of previous _N ^{p n-1} and the calculation area volume _{V d (= (2πL 0)} 3), collision in n-1/2 step frequency _{N c} ^{n-1 / 2} is calculated as follows.

This collision frequency and n-1/2 particle number in the step density ^{_{n p n-1/2 (}} = 0.5 (N p n-1 + N p n) / V d), collision in n-1/2 step The frequency factor is calculated by the following formula.

The average collision factor is calculated by averaging the collision frequency factor over time.

If there are N _p particles in the calculation area, N _p (N _p −1) / ^{2 to} 0.5 N _p ^two collision determination calculations must be performed unless a special method is used. Cell index method (Allen, MP, and Tildesley, DJ, “Computer Simulation of Liquids”, Oxford University Press, (1987), p. 408) used in molecular dynamics as a method to reduce the number of collision determinations . As shown in FIG. 3, the collision frequency region is subdivided into cells, and a cell including a certain particle and particles in a total of 27 cells before and after the x, y, and z directions surrounding the cell. Determine collision. However, in the case of sequentially determining collisions, instead of 27 cells, 8 cells (for example, cells including target particles and cells ahead in the x, y, and z directions) are determined. If this is the case, there is no duplication judgment and efficiency is improved. In this cell index method, a number is assigned to each particle, and a list (cell index) for managing what number and what number particle is in which cell is created.

As shown in FIG. 3, when the calculation region is divided into ^three _cells of M _cell , the average number of particles N _{p, cell} in each _cell is N _{p, cell} = N _p / M _cell ³ . The number of collision determinations for one particle is the number of collision determinations (Np _{, cell} −1) / 2 with other particles in the cell containing the particles, and the in-cell particles in the x, y, z direction rear side. The sum of the collision determination times 7N _{p, cell} , that is, (N _{p, cell} −1) / 2 + 7N _{p, cell} is estimated. After all, the total number of collision determination count is a _{N p × ((N p,} cell -1) / 2 + 7N p, cell) ~7.5N p 2 / M cell 3. That is, if M _cell > 2.5, the number of collision determinations can be reduced. In this estimation, the larger the M _cell , the smaller the number of collision determinations. However, in a situation where there are almost no particles in the cell, the process of checking whether there is a particle in the cell is wasted. Therefore, it is desirable to set M _cell so that N _{p, cell} is 1 or more.

In the present embodiment, gas phase turbulence is calculated using a finite difference method, and collision between particles moving in the turbulent flow is detected. In order to support distributed memory type parallel computation using the MPI library, the following three corresponding MPI communication functions must be realized. That is, an MPI communication program corresponding to the function must be described.
(A) Acquisition of gas phase data of sleeve region necessary for difference calculation (B) Processing to move particles protruding from MPI division region (calculation region assigned to each node 10) to the next process (C) MPI When a particle pair sharing particle information in the vicinity of the boundary of the divided region is close to the boundary of the MPI divided region, the collision may not be detected unless the processes on both sides share the information. In order to deal with this problem, the above operation (C) is required. The present invention is for efficiently processing (B) and (C) by utilizing a cell index.

  Subsequently, the function according to the present invention of the node 10 will be described in detail. As shown in FIG. 1, the node 10 includes a particle state calculation unit 11, a period determination unit 12, an own virtual particle information output unit 13, another virtual particle information input unit 14, an own virtual particle state output unit 15, The other virtual particle state input unit 16 and the output unit 17 are provided. Note that more specific functions of each component and more specific processing performed by each component will be described in the description of the processing flow described later. These functions are realized, for example, by executing a program in the node 10.

  The particle state calculation unit 11 is a particle state calculation unit that calculates the state of particles in its own region at each time step. The own area is a calculation area assigned to the own node 10. Specifically, the particle state calculation unit 11 calculates the particle state by the method described above. Parameters necessary for the calculation are input from a device that comprehensively controls the above-described distributed memory type parallel calculation. Further, the calculation formulas and the like necessary for the calculation are stored in advance in the particle state calculation unit 11 and are read and calculated (the same applies to other components).

  The particle state calculation unit 11 performs processing according to the cell index method as described above for the determination and calculation of the interaction between particles (in this embodiment, collision between particles). That is, each particle is given a number that is information for identifying the particle. The particle state calculation unit 11 stores and manages the number of the particle in its own area in association with the information specifying the cell containing the particle. As described above, the cell is a region obtained by further dividing the self region.

  As will be described later, the particle state calculation unit 11 specifies information on other region virtual particles input by the other virtual particle information input unit 14 and the state of the other region virtual particles input by the other virtual particle state input unit 16. The state of the particles is calculated using information indicating A virtual particle is a particle whose state is to be calculated in its own region (own node 10), but is used to calculate the state of a particle in an adjacent region (another node 10) that is a divided calculation region adjacent to its own region. Necessary particles. The virtual particle information is transmitted and received between the nodes 10 and the information is used for the calculation of the particles. Among the particles in the self region, virtual particles used for calculation in the adjacent region are referred to as self region virtual particles. Among particles in an adjacent region, virtual particles used for calculation in the own region are called other region virtual particles.

  The particle state calculation unit 11 also transmits / receives information on particles that protrude from the own region and enter the adjacent region and information on particles that protrude from the adjacent region and enter the own region with the adjacent region (another node 10). Do.

  Further, as described above, the particle state calculation unit 11 calculates the state of the fluid by the finite difference method, and calculates the state of the particles in consideration of the influence of the fluid.

  The particle state calculation unit 11 specifies the proximity relationship between the particles based on the information specifying the other region virtual particles input by the other virtual particle information input unit 14, and determines the particle state based on the proximity relationship. calculate. In this embodiment, the proximity relationship is a particle pair for which a collision is determined. Specifically, a cell index is used as will be described later.

  The period determination unit 12 is a period determination unit that determines whether the time step of the calculation in the particle state calculation unit 11 is a time step of a preset period. A value indicating this period is a natural number of 2 or more, and is stored in the period determination unit 12 in advance. The period determination unit 12 makes the above determination by referring to information on the current time step managed by the particle state calculation unit 11, for example. When the period determination unit 12 determines that the time step of the calculation is a time step of a preset period, the period determination unit 12 notifies the particle state calculation unit 11 and the own virtual particle information output unit 13 to that effect.

  The own virtual particle information output unit 13 specifies the own region virtual particle at the time step determined to be the time step of the cycle set in advance by the cycle determination unit 12, and indicates the identified own region virtual particle Is a self-virtual particle information output means for outputting to the other node 10 for calculating the state of particles in the adjacent region. The identification of the self-region virtual particles is performed based on the position of the particles calculated by the particle state calculation unit 11. Of the cells in the own region, the particles located in the cells adjacent to the adjacent region are set as the own region virtual particles. The own virtual particle information output unit 13 specifies the own region virtual particle for each adjacent region, and outputs the number of the own region virtual particle to another node 10 that calculates the state of the particles in the adjacent region. The own virtual particle information output unit 13 notifies the own virtual particle state output unit 15 of the number of the own region virtual particle.

  The other virtual particle information input unit 14 is another virtual particle information input means for inputting the number of the other region virtual particle that is information for specifying the other region virtual particle from another node 10 that calculates the state of the particle in the adjacent region. is there. The other virtual particle information input unit 14 notifies the particle state calculation unit 11 of the input number of the other region virtual particles.

  The own virtual particle state output unit 15 is a self virtual particle state output unit that outputs information indicating the state of the own region virtual particle calculated by the particle state calculation unit 11 at each time step to another node 10. The own virtual particle state output unit 15 obtains the particle state information indicated by the number of the own region virtual particle notified from the own virtual particle information output unit 13 from the particle state calculation unit 11 and outputs the information to another node 10. To do. The node 10 to be output is the node 10 to which the virtual particle information output unit 13 outputs the number information. The information on the state of the particles to be output includes information indicating the position (coordinates) of the particles and information indicating the velocity (direction and size) of the particles.

  The other virtual particle state input unit 16 is another virtual particle state input unit that inputs information indicating the state of the other region virtual particles from another node 10 at each time step. The own virtual particle state output unit 15 outputs information indicating the state of the input other region virtual particles to the particle state calculation unit 11. The other region virtual particles to which the particle information is input are other region virtual particles to which information specified by the other virtual particle information input unit 14 is input. The other virtual particle state input unit 16 outputs information indicating the state of the input other region virtual particles to the particle state calculation unit 11.

  The output unit 17 is an output unit that outputs information according to the calculation by the particle state calculation unit 11. Specifically, for example, the output unit 17 outputs the information for calculating the above-described average collision frequency factor obtained by the calculation by the particle state calculation unit 11 to an apparatus that centrally controls the distributed memory type parallel calculation. To do.

  Subsequently, a process (particle state calculation method) executed by the node 10 according to the present embodiment will be described using the flowcharts of FIGS. 4 and 5. First, an outline of processing of the entire particle state calculation system 1 will be described using the flowchart of FIG. 4, and then processing according to the present invention by each node 10 will be described using the flowchart of FIG.

  In the particle state calculation system 1, first, parameters for calculation are input to an apparatus that comprehensively controls distributed memory type parallel calculation (S01). This parameter includes the size of the calculation area, the number of grids for fluid calculation, the number of cells for particle calculation, the number of particles, the size of each particle, the density, the speed, the number of divisions for distributed memory type parallel calculation, and the time increment. This is information such as (time step interval). This input is performed, for example, by an operation of the device by the user.

  In the apparatus for comprehensively controlling the distributed memory type parallel calculation, the calculation area assigned to each node 10 is determined based on the input parameters, and information necessary for the calculation is output to each node (S02). ). Each node 10 performs calculation (S03). Subsequently, information (for example, average collision frequency factor) corresponding to the calculation result of each node 10 is output (S04). As for the output of the results, for example, the calculation results by each node 10 are collected in a device that comprehensively controls the distributed memory type parallel calculation, and the final result is generated and output by the device.

Subsequently, processing in each node 10 (corresponding to S03 in the flowchart of FIG. 4) will be described using the flowchart of FIG. FIG. 6 shows calculation areas corresponding to the respective processes. In the diagram shown in FIG. 6, the calculation area is a two-dimensional area, and M _{cell, x} × M _{cell, y} cells (M _{cell, x} = M _cell / M _x and M _{cell, y} = M) that a certain process takes charge of. _cell / M _y ) is drawn with a solid line. A region drawn with a solid line is also an MPI division region (own region) that the target process is in charge of, and a dotted line portion indicates a region in which another process (neighboring process) is in charge of calculation. Hereinafter, the particles in the MPI divided region at the time of creating the cell index are referred to as real particles. For collision determination, it is necessary to consider particles in the region of M _{cell, x} + 1 and M _{cell, y} + 1. Those particles are virtual particles. Virtual particles are used only for collision detection and their motion is not calculated. In FIG. 6, the black dots are real particles, and the white dots are (other) virtual particles.

  In this process, first, the particle state calculation unit 11 calculates the fluid state at a predetermined time step (S11, particle state calculation step). The calculated fluid state is the fluid velocity and pressure. In addition, the particle state calculation unit 11 calculates the state of the actual particle in a predetermined time step (S12, particle state calculation step). The calculated real particle state is the position and velocity of the real particle. These calculations are performed by the Runge-Kutta method, as described above, based on the fluid and actual particle states in the previous time step or initial state and the fluid state after the calculation. It should be noted that the calculation of the state of the fluid and the actual particles is not necessarily performed sequentially, and is calculated alternately or simultaneously for each of the fluid and the actual particles, such as the primary step and the secondary step of the Runge-Kutta method. The state shown in FIG. 6A is the previous time step or initial state (that is, immediately before S11), and the state shown in FIG. 6B is the state after calculation (that is, immediately after S12).

  After S12, as shown in FIG. 6B, (real) particles appear from the MPI division calculation area. If the processing of the protruding particles (the processing of (B) described above) is executed step by hour, a large amount of communication between the nodes 10 may occur, which may reduce parallel performance. Therefore, processing of particles that protrude only in a preset period is performed as follows. In order to continue the calculation of the protruding actual particles, fluid velocity data outside the MPI division calculation area is required. Each process (particle state calculation unit 11) has sleeve data for four grids for the fourth-order center difference calculation in accordance with the processing of (A) described above. Therefore, u is calculated without extra MPI communication (communication between nodes 10) as long as the particles do not protrude beyond 3.5Δ (3.5Δ instead of 4Δ because a staggered lattice is used). can do. In other words, it is an efficient technique that is possible because the fluid calculation is performed by the finite difference method. In addition, as shown below, cell index creation is also performed at preset intervals.

  After S12, the cycle determination unit 12 determines whether or not the current time step is a preset time step (S13, cycle determination step). When it is determined that the current time step is a time step of a preset cycle (Yes in S13), the cycle determination unit 12 notifies the particle state calculation unit 11 and the own virtual particle information output unit 13 to that effect. In this case, the processing (S14 to S19) is performed.

  First, as shown in FIG.6 (c), the information regarding the other virtual particle memorize | stored in the particle | grain state calculation part 11 is removed (S14). Subsequently, the particle state calculation unit 11 inputs / outputs the information of the actual particles that have exited from the own region (S15, particle state calculation step). Specifically, the information on the real particles that have exited from the own region is output to a node that is responsible for calculating the region in which the actual particle is located, and is deleted from the actual particles in the own region. This real particle becomes a real particle in another region. On the other hand, the information of the particles that have exited from the other area and entered the own area is received from the node 10 that is responsible for the calculation of the other area, and becomes the actual particles in the own area (particles indicated by arrows in FIG. Is). That is, delivery between nodes (regions) of real particles is performed. Here, communication (MPI communication) between the nodes 10 occurs.

  Subsequently, a cell index is generated (updated) by the particle state calculation unit 11. The cell index is information in which the number of an actual particle included in the cell (located on the cell) is associated with information for specifying the cell (S16, particle state calculation step). The list of real particles in the cell (2, 2) shown in FIG. 6 (e) is [4, 102, 605, 651]. The cell index is stored and managed by the particle state calculation unit 11. The cell (2, 2) indicates the second cell on the x-axis and the second cell on the y-axis (the same applies hereinafter).

  Subsequently, the self-virtual particle information output unit 13 identifies the self-region virtual particle and outputs information indicating the self-region virtual particle to another node 10 that calculates the state of the particle in the adjacent region (S17, self Virtual particle information output step). Specifically, in the cell index by the particle state calculation unit 11, real particles associated with a cell in which the cell included in the cell is determined to be the own region virtual particle are set as the own region virtual particle. Such a cell is determined by the positional relationship between the regions. Information identifying the cell and the number of the own region virtual particle are output from the own virtual particle information output unit 13 to another node 10. The self-virtual particle number is notified from the self-virtual particle information output unit 13 to the self-virtual particle state output unit 15.

  On the other hand, information in which the number of the other region virtual particle transmitted from another node 10 and the cell including the particle are associated is input by the other virtual particle information input unit 14 (S18, other virtual particle). Information input step). This information is transmitted by a process corresponding to S17 in another node 10. The above information is input from the other virtual particle information input unit 14 to the particle state calculation unit 11. Thereby, as shown in FIG.6 (f), the number of the other virtual particle located in the cell adjacent to an own area | region is grasped | ascertained by the particle | grain state calculation part 11. FIG. In S17 and S18, communication (MPI communication) between the nodes 10 occurs.

  Subsequently, the particle state calculation unit 11 generates a particle pair for checking the presence / absence of collision for each particle in its own region (S19, particle state calculation step). This particle pair is a pair of a target particle, a cell including the target particle, and a particle located in a cell adjacent to the cell. In generating this pair, other virtual particles are also considered (a pair of real particles and other virtual particles).

  The (i, j) -th cell ((i, j) cell) includes itself (i, j), (i, j ± 1), (i ± 1, j), (i ± 1). , J ± 1) are surrounded by a total of nine cells (27 in the case of three dimensions). However, when checking the collision by sequentially increasing the cell number, a total of four (i, j), (i, j + 1), (i + 1, j), (i + 1, j + 1) (8 in the case of three dimensions). ) Cells, you can check all combinations without duplication. Therefore, when creating a list of particle pairs, it is not necessary to include cells with i = 0 or j = 0. The list in which the number 130 is a pair shown in FIG. 6 (f) is [332, 69, 871, 777, 109, 4, 605, 102, 651] which are particles contained in the above four cells.

  After S19 or after S12, when the cycle determination unit 12 determines that the current time step is not a preset time step (No in S13), the following processing is performed.

  First, the self-virtual particle state output unit 15 acquires information indicating the state of the current time step of the self-region virtual particles notified from the self-virtual particle information output unit 13. The acquired information is the information calculated by the particle state calculation unit 11 in S12. The acquired information is output from the own virtual particle state output unit 15 to the other node 10 (S20, own virtual particle state output step).

  On the other hand, information indicating the state of other region virtual particles transmitted from another node 10 is input by the other virtual particle state input unit 16 (S21, other virtual particle state input step). This information is transmitted by a process corresponding to S20 in another node 10. The above information is input from the other virtual particle state input unit 16 to the particle state calculation unit 11. Thereby, as shown in FIG. 6G, the particle state calculation unit 11 grasps the position and velocity that are the states of the other virtual particles located in the cell adjacent to the own region. In S21 and S22, communication (MPI communication) between the nodes 10 occurs.

  Subsequently, as shown in FIG. 6H, the particle state calculation unit 11 determines the collision between particles as described above, and performs the calculation based on the determination (S22, particle state calculation step). ). This determination is made among all the pairs generated in S19. Judgment of collision between real particles is performed based on the state of particles calculated in S12. The determination of the collision between the real particle and the other region virtual particle is performed based on the particle state calculated in S12 for the real particle and the particle state input in S21 for the other region virtual particle.

  Subsequently, the particle state calculation unit 11 determines whether to end the calculation (S23, particle state calculation step). This determination is made based on a determination criterion set or input in advance in the particle state calculation unit 11. When it is determined that the calculation is not completed (No in S23), the process proceeds to the next time step and the processing from S11 is repeated (S24, particle state calculation step).

  When it is determined that the calculation is finished (Yes in S23), the output unit 17 outputs information according to the calculation by the particle state calculation unit 11 (S25, output step). Thus, the process in the node 10 ends. In addition, after updating the time step in the above processing (S24), it is also possible to calculate and output statistics such as particle collision, particle velocity and particle energy once every step or several steps. Good.

In the above-described processing, once a cell index is generated, it is not updated until a time step of a preset period elapses, so that the particles may move from the original cell during that time. However, if the moving distance of the particles during the time step of the above period is smaller than the cell width Δ _cell (= 2πL ₀ / M _cell ), the possibility of leaking a particle pair that may collide can be ignored. In the present embodiment, for example, when Δ _cell = 4Δ (that is, M _cell = N / 4), when the above period is set to 8 and the Courant number is set to 0.2 or less, the particles are about the fluid velocity. Considering movement, the particles move at most about 1.6Δ during the time step of the above period. In other words, in this setting, collision detection leaks are negligibly small. In fact, in some cases, it was confirmed that there was no change in the collision pair when the period was set to 1 and when it was set to 8.

  As described above, in this embodiment, the particles related to information exchanged between the nodes 10 are limited virtual particles, so that necessary information can be surely input in each node 10 and communication between the nodes 10 can be efficiently performed. To do. Further, communication between the nodes 10 becomes more efficient by exchanging information for specifying virtual particles as a time step having a preset period. As a result, according to this embodiment of the present invention, it is possible to efficiently calculate the particle state by distributed memory type parallel calculation, and to calculate a larger model than the result.

  Moreover, if the particles in the fluid are calculated as in the present embodiment, a large-scale model calculation related to the particles in the fluid can be performed. However, it is not always necessary to calculate the particles in the fluid, and can be applied to those that calculate the particle state by distributed memory type parallel calculation.

  In addition, as in this embodiment, if the cell pair is used to determine the interaction between pairs of particles, it is possible to reliably calculate a model in which interactions between particles such as collisions between particles occur. Can do.

In addition, about calculation of the particle | grains in embodiment mentioned above, it is good also as mixing a particle | grain after a steady gaseous-phase turbulent flow field is formed. In addition, the collision may be determined after a certain period of time has elapsed since the mixing of particles (for example, T ₀ (= L ₀ / U ₀ ) or more).

  In the present embodiment, an example of simulating a cloud has been shown. However, the present invention is not necessarily limited to calculation for a cloud, and can be applied as long as it calculates a particle state. For example, it can be applied to asteroid generation simulation.

The specific calculation result by embodiment mentioned above is shown below. In addition, after confirming that a steady gas phase turbulent phase field was formed, particles were mixed. When calculating particle motion, the time step Δt was set to 0.5τ _p or less. Δt that satisfies both this condition and the condition based on the Courant number of the fluid calculation was used for both the fluid calculation and the particle calculation. Further, the collision determination was started after T ₀ (= L ₀ / U ₀ ) or more elapsed from the mixing of particles.

The following table shows the particle calculation conditions. The case names shown in the leftmost column correspond to the above table. In all cases, the particle number density _np was set to 10 ⁸ order, which is a typical value of clouds. The average number of particles contained in the cell was 1 for N64 and N256, and 4.8 for N512.

The following table shows the turbulent characteristic quantities of the gas-phase isotropic turbulent flow field that statistically reached a steady state. In the table, u ′ is the velocity fluctuation intensity, Re _λ (= u′l _λ / ν, l _λ is Taylor microscale) is Taylor microscale reference turbulent Reynolds number, and k _max (= N / 2) is the maximum effective wave number. , L _η (= (ν ³ / ε) ^1/4 , ε is the energy discrete rate) is the Kolmogorov scale. k _max l _η is an index of calculation resolution, and is said to be about 1 when discussing two-dimensional quantities such as kinetic energy and about 2 when discussing high dimensions such as skewness and flatness. Yes. In the conventional method, the particle collision calculation is often set to about 1.3 to 1.4, but in this example, the resolution is set to about 2, which is higher than usual.

FIG. 7 shows the result of the collision frequency factor in the case of N64 (Re _λ = 52.3). The vertical axis is the collision frequency factor made dimensionless by the local velocity gradient λ (= (ε / ν) ^1/2 ) and the collision radius R (= 2r), and the horizontal axis is the particle radius made dimensionless on the Kolmogorov scale. is there. The average collision frequency factor in each run was calculated with the number of particles reduced by 5% due to collision as one run. For comparison, Onishi et al. (2009) (Onishi, R., Takahashi, K. and Komori, S., “Influence of Gravity on Collisions of Monodispersed Droplets in Homogeneous Isotropic Turbulence”, Physics of Fluids, Vol. 21 ( 2009) and 125108), the results for Re _λ = 54.3 are also shown. Onishiet al. (2009) used a pseudospectral method for fluid calculations. This is different from the present embodiment using the finite difference method, but the setting of the size of the calculation area, the physical property value, the number of calculation grids, the number of particles, etc., and the calculation method of the standard deviation are the same as in this embodiment. Since the result of this embodiment and the result of Onishi et al. (2009) agree well, the reliability of this embodiment was confirmed.

Next, the collision frequency factor will be described. FIG. 8 shows the collision frequency factor of a particle having a Stokes number St (= τ _p λ) of 1. The vertical axis represents the dimensionless collision frequency factor, and the horizontal axis represents the turbulent Reynolds number Re _λ . The black plot represents the DNS results. Circle plots are according to this embodiment, triangle plots are according to Onishi et al. (2009), and square plots are according to Wottitez et al. (2009) (non-patent document 1). White plots are shown by Onishi et al. (2007) (Ryo Onishi, Satoru Komori, “Modeling of collision factors between particles of the same diameter in turbulent flow”, Journal of the Japan Society of Mechanical Engineers, B, Vol. 73, No. 2007), pp. 1307-1314.). The arrow extending from the model indicates the value that the model is trying to predict.

In the case of a low Reynolds number of Re _λ <85, which has been obtained in previous studies, the model has a few percent error, but it can be well predicted that the collision frequency factor tends to decrease as the Reynolds number increases. However, in the case of Re _λ > 85 high Reynolds number obtained for the first time by this embodiment, the model prediction error is large. While the model predicts that the collision frequency factor decreases rapidly with increasing Reynolds number, the DNS results do not show such a clear decreasing trend. This is because the model of Onishi et al. (2007) uses empirical parameters obtained based on data in the case of a low Reynolds number.

  For example, Falkovich et al. (2002) (Falkovich, G., Fouxon, A. and Stepanov, MG, “Acceleration of RainInitiation by Cloud Turbulence”, Natrure, Vol.419 (2002), pp.151-154) It is predicted that the “sling effect” caused by turbulence intermittency will increase the collision frequency when the Reynolds number is high. Since such effects are not taken into account, Onishi et al. (2007) model underestimates the collision frequency for high Reynolds numbers.

  Next, parallel computing performance will be described. Large scale parallel computation requires high linear scalability. For the evaluation of linear expandability, two types of strong linear expandability (strong scaling) and weak linear expandability (weak scaling) are used. In order to evaluate the former, after fixing the entire calculation scale, it is examined how much the calculation time is reduced when the number of calculation processes (the number of nodes 10 in the present embodiment) is increased. For example, when the number of calculation processes is doubled, if the calculation time is halved, it can be said that it has ideal strong linear expandability. On the other hand, in order to evaluate the latter, the calculation scale for the calculation process is fixed, and how much the calculation time changes when the number of calculation processes is increased is examined. For example, even if the overall calculation scale and the number of calculation processes are both doubled, it can be said that the calculation has a weak weak linear extensibility if the calculation time does not change.

FIG. 9 shows (a) the result of the calculation time [s] required for one time step and (b) the result of the calculation speed [MFLOPS] when the number of calculation processes is changed with respect to the N256 calculation in this embodiment. Indicates. The horizontal axis in both figures is the number of arithmetic processes. This is a result of examining the above-described strong linear expandability because the number of operation processes was changed after fixing the entire calculation scale (256 ³ lattices, 262, 144 particles). The SGI ALTIX 4700 system (6.4 GFLOPS / core) was used for this performance measurement. From FIG. 9A, it can be seen that the computation time decreases in inverse proportion to the number of computation processes. Correspondingly, in FIG. 9B, it can be seen that the calculation speed increases almost in proportion to the number of processes.

  For example, a large-scale parallel DNS for an isotropic turbulent flow field using the IBM BlueGene system has been reported to have a performance of about 10% of the theoretical peak performance. In the present embodiment, high performance equivalent to that of fluid calculation alone can be achieved even after calculation of particle motion and collision determination.

As described above, the present embodiment is a high-efficiency large-scale parallel calculation method that acquires a collision data between particles by calculating a turbulent flow field by a finite difference method, calculating a particle motion by a Lagrangian tracking method. This embodiment has the following features.
(I) As a compulsory method for obtaining a steady isotropic turbulent flow field, an RCF method that does not impair the high parallelization rate of the finite difference method was used.
(Ii) The cell index method was used to efficiently detect the collision between particles.
(Iii) Both shared memory type parallel (automatic parallelization library) and distributed memory type parallel (MPI library) are supported, and hybrid calculation using both simultaneously is supported.
(Iv) The code is designed so that the communication of the particle information and the update of the cell index information protruding from the MPI parallel calculation area is performed only at each time step of a predetermined cycle, and the inter-process communication cost in the distributed memory type parallel calculation is reduced. Reduced.

Actually, by executing the processing based on this embodiment, collision frequency data of inertial particles at a high Reynolds number of Re _λ = 209 could be obtained. This is valuable data that greatly updates the maximum Reynolds number Re _λ = 85 in the conventional particle collision calculation. In addition, the existing collision frequency model was also verified using the collision frequency data at the high Reynolds number obtained in this embodiment. As a result, it became clear that the existing collision frequency model underestimates the collision frequency when the Reynolds number is high.

DESCRIPTION OF SYMBOLS 1 ... Particle state calculation system, 10 ... Node, 11 ... Particle state calculation part, 12 ... Period determination part, 13 ... Self virtual particle information output part, 14 ... Other virtual particle information input part, 15 ... Self virtual particle state output part 16, other virtual particle state input unit, 17 ... output unit.

Claims (4)

A particle state calculation device that repeatedly calculates a state of particles arranged in a self region that is one of a plurality of divided calculation regions for each time step,
Particle state calculating means for calculating the state of the particles in the local region at each time step;
A period determining means for determining whether the target time step in the particle state calculating means is a time step of a preset period;
Based on the position of the particle calculated by the particle state calculating means, the particle adjacent to the self-region based on the position of the particle calculated by the particle state calculating means at the time step determined to be a time step of a period preset by the period determining means. The self-region virtual particle used for the calculation in the adjacent region that is the divided calculation region is specified, and the information indicating the specified self-region virtual particle is supplied to another particle state calculation device that calculates the state of the particle in the adjacent region. Self-virtual particle information output means for outputting;
Other virtual particle information input means for inputting information specifying other region virtual particles used for calculation in the own region from the other particle state calculation device;
A self-virtual particle state output unit that outputs information indicating the state of the self-region virtual particle calculated by the particle state calculation unit in each time step to the other particle state calculation device;
Other virtual particle state input means for inputting information indicating the state of the other region virtual particles from the other particle state calculation device at each time step;
Output means for outputting information according to the calculation by the particle state calculation means,
The particle state calculation means includes information specifying the other region virtual particles input by the other virtual particle information input means, and information indicating the state of the other region virtual particles input by the other virtual particle state input means. A particle state calculation device that calculates the state of particles using a particle. The particles are particles in a fluid;
The particle state calculation means calculates the state of the fluid by a finite difference method, and calculates the state of the particles in consideration of the influence of the fluid.
The particle state calculation apparatus according to claim 1.   The particle state calculation means specifies the proximity relationship between particles based on the information specifying the other region virtual particles input by the other virtual particle information input means, and the particle state based on the proximity relationship The particle state calculation apparatus according to claim 1 or 2, wherein A particle state calculation method by a particle state calculation device that repeatedly calculates a state of particles arranged in its own region which is one of a plurality of divided calculation regions for each time step,
A particle state calculating step for calculating the state of the particles in the local region at each time step;
A period determination step for determining whether or not the target time step in the particle state calculation step is a time step of a preset period;
Based on the position of the particle calculated in the particle state calculation step, the particle adjacent to the self-region based on the position of the particle calculated in the particle state calculation step is the time step determined to be a time step of a period set in advance in the cycle determination step. The self-region virtual particle used for the calculation in the adjacent region that is the divided calculation region is specified, and the information indicating the specified self-region virtual particle is supplied to another particle state calculation device that calculates the state of the particle in the adjacent region Self virtual particle information output step to output;
Other virtual particle information input step for inputting information specifying other region virtual particles used for calculation in the own region from the other particle state calculation device;
A self-virtual particle state output step of outputting information indicating the state of the self-region virtual particle calculated in the particle state calculation step in each time step to the other particle state calculation device;
Other virtual particle state input step of inputting information indicating the state of the other region virtual particles from the other particle state calculation device in each time step;
An output step for outputting information according to the calculation in the particle state calculation step,
In the particle state calculation step, information specifying the other region virtual particles input in the other virtual particle information input step, and information indicating the state of the other region virtual particles input in the other virtual particle state input step A particle state calculation method for calculating the state of particles using a particle.