PeTar
N-body code for collisional gravitational systems
cluster_list.hpp
Go to the documentation of this file.
1 #pragma once
2 #include<particle_simulator.hpp>
3 #include<unordered_map>
4 #include<map>
5 #include"ptcl.hpp"
6 #include"Common/binary_tree.h"
7 
8 //extern const PS::F64 SAFTY_OFFSET_FOR_SEARCH;
9 
10 #ifndef ARRAY_ALLOW_LIMIT
11 #define ARRAY_ALLOW_LIMIT 1000000000
12 #endif
13 
14 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
15 template <class Ttree>
16 void SetRankComm(const PS::F64ort pos_domain[], Ttree &tree,
17  PS::ReallocatableArray<PS::S32> & rank_neighbor,
18  const PS::F64 eps_sq = 0.0){
19  rank_neighbor.clearSize();
20 // static const PS::F64 SAFTY_FACTOR_FOR_SEARCH_SQ = 1.01;
21  const PS::S32 my_rank = PS::Comm::getRank();
22  const PS::S32 n_proc = PS::Comm::getNumberOfProc();
23 // const PS::F64ort pos_my_domain = pos_domain[my_rank];
24  PS::F64ort pos_my_domain = tree.getOuterBoundaryOfLocalTree();
25  PS::F64ort pos_my_domain_in = tree.getInnerBoundaryOfLocalTree();
26  PS::F64ort * pos_domain_out = new PS::F64ort[n_proc];
27  PS::Comm::allGather(&pos_my_domain, 1, pos_domain_out);
28  //MPI_Allgather(&pos_my_domain, 1, PS::GetDataType<PS::F64ort>(),
29  // pos_domain_out, 1, PS::GetDataType<PS::F64ort>(), MPI_COMM_WORLD);
30  PS::F64ort * pos_domain_in = new PS::F64ort[n_proc];
31  PS::Comm::allGather(&pos_my_domain_in, 1, pos_domain_in);
32  //MPI_Allgather(&pos_my_domain_in, 1, PS::GetDataType<PS::F64ort>(),
33  // pos_domain_in, 1, PS::GetDataType<PS::F64ort>(), MPI_COMM_WORLD);
34 // const PS::F64 r_crit_sq = (EPJSoft::r_search * EPJSoft::r_search + eps_sq)*SAFTY_FACTOR_FOR_SEARCH_SQ*1.01;
35 // pos_domain[my_rank].getOuterBoundaryOfLocalTree
36  rank_neighbor.clearSize();
37  for(int i=0; i<n_proc; i++){
38  if(i==my_rank) continue;
39  //else if( r_crit_sq >= pos_my_domain.getDistanceMinSQ(pos_domain[i]) ){
40  else if( pos_my_domain.contained(pos_domain_in[i]) || pos_my_domain_in.contained(pos_domain_out[i])){
41  rank_neighbor.push_back(i);
42  }
43 // rank_neighbor.push_back(i);
44  }
45  delete [] pos_domain_in;
46  delete [] pos_domain_out;
47 }
48 #endif
49 
50 struct Cluster{
51  PS::S32 id_;
52  PS::S32 n_ptcl_; // include outer particle
53  PS::S32 n_ptcl_stored_; // exclude outer particle
54  PS::S32 adr_head_; // for cluster_isolated, adr_sys_multi_cluster_isolated_
55  PS::S32 rank_;
56  Cluster(): id_(-1), n_ptcl_(0), n_ptcl_stored_(0), adr_head_(-1), rank_(-1){}
57  Cluster(const PS::S32 _id, const PS::S32 _n_ptcl, const PS::S32 _n_ptcl_stored,
58  const PS::S32 _adr_head, const PS::S32 _rank):
59  id_(_id), n_ptcl_(_n_ptcl), n_ptcl_stored_(_n_ptcl_stored),
60  adr_head_(_adr_head), rank_(_rank){}
61  void clear(){
62  id_ = -1; n_ptcl_ = 0; n_ptcl_stored_ = 0; adr_head_ = -1; rank_ = -1;
63  }
64  void dump(){
65  std::cout<<" id="<<id_<<" n_ptcl="<<n_ptcl_
66  <<" n_ptcl_stored="<<n_ptcl_stored_
67  <<" adr_head="<<adr_head_
68  <<" rank="<<rank_<<std::endl;
69  }
70 };
71 
72 struct PtclCluster{
73  PS::S32 id_;
74  PS::S32 adr_sys_;
75  PS::S32 adr_ngb_head_; // adr of id_ngb_multi_cluster
76  PS::S32 n_ngb_;
77  bool flag_searched_;
78  PtclCluster * next_;
79  PS::S32 rank_org_;
80  PS::S32 id_cluster_;
81  PtclCluster(): id_(-10), adr_sys_(-1), adr_ngb_head_(-1), n_ngb_(0),
82  flag_searched_(false), next_(NULL), rank_org_(-1), id_cluster_(id_){}
83  PtclCluster(const PS::S32 _id, const PS::S32 _adr_sys, const PS::S32 _adr_ngb_head,
84  const PS::S32 _n_ngb, const bool _flag_searched, PtclCluster * const _next,
85  PS::S32 _rank_org):
86  id_(_id), adr_sys_(_adr_sys), adr_ngb_head_(_adr_ngb_head),
87  n_ngb_(_n_ngb), flag_searched_(_flag_searched), next_(_next), rank_org_(_rank_org),
88  id_cluster_(_id){}
89  void clear(){
90  id_ = -10; adr_sys_ = -1; adr_ngb_head_ = -1; n_ngb_ = 0;
91  flag_searched_ = false; next_ = NULL; rank_org_ = -1; id_cluster_ = id_;
92  }
93  void dump(){
94  std::cout<<" id="<<id_<<" adr_sys="<<adr_sys_<<" adr_ngb_head_="<<adr_ngb_head_
95  <<" n_ngb="<<n_ngb_<<" flag_searched_="<<flag_searched_
96  <<" rank_org="<<rank_org_<<" id_cluster_="<<id_cluster_<<std::endl;
97  }
98  void setIdClusterImpl(PS::ReallocatableArray<PtclCluster> & ptcl,
99  PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > & adr_ngb){
100  for(PS::S32 i=0; i<n_ngb_; i++){
101  id_cluster_ = (ptcl[adr_ngb[adr_ngb_head_+i].second].id_cluster_ < id_cluster_)
102  ? ptcl[adr_ngb[adr_ngb_head_+i].second].id_cluster_ : id_cluster_;
103  }
104  }
105 };
106 
107 class PtclOuter{
108 public:
109  PS::S32 id_;
110  PS::S32 id_ngb_;
111  PS::S32 rank_org_;
112  PtclOuter(): id_(-10), id_ngb_(-1), rank_org_(-1){}
113  PtclOuter(const PS::S32 _id, const PS::S32 _id_ngb, const PS::S32 _rank_org):
114  id_(_id), id_ngb_(_id_ngb), rank_org_(_rank_org){}
115  void clear(){
116  id_ = -10; id_ngb_ = -1; rank_org_ = -1;
117  }
118  void dump(){
119  std::cout<<" id="<<id_<<" id_ngb_="<<id_ngb_<<" rank_org_="<<rank_org_<<std::endl;
120  }
121 };
122 
123 class Mediator{
124 public:
125  PS::S32 id_;
126  PS::S32 adr_sys_;
127  PS::S32 adr_pcluster_;
128  PS::S32 id_cluster_;
129  PS::S32 rank_send_;
130  Mediator():id_(-10), adr_sys_(-1), adr_pcluster_(-1), id_cluster_(-1), rank_send_(-1){}
131  Mediator(const PS::S32 _id, const PS::S32 _adr_sys,
132  const PS::S32 _adr_pcluster, const PS::S32 _id_cluster,
133  const PS::S32 _rank_send):
134  id_(_id), adr_sys_(_adr_sys), adr_pcluster_(_adr_pcluster), id_cluster_(_id_cluster),
135  rank_send_(_rank_send){}
136  void clear(){
137  id_ = -10; adr_sys_ = -1; adr_pcluster_ = -1; id_cluster_ = -1; rank_send_ = -1;
138  }
139  void dump(){
140  std::cout<<" id="<<id_<<" adr_sys="<<adr_sys_
141  <<" adr_pcluster_="<<adr_pcluster_<<" id_cluster="<<id_cluster_
142  <<" rank_send_="<<rank_send_<<std::endl;
143  }
144 };
145 
146 class PtclComm: public Ptcl{
147 public:
148  PS::S32 id_cluster;
149 
150  template <class Tp>
151  PtclComm(const Tp &p): Ptcl(p), id_cluster(-1) {}
152  PtclComm() {}
153 
154  void print(std::ostream & fout){
155  Ptcl::print(fout);
156  fout<<" id_cluster="<<id_cluster<<std::endl;
157  }
158 };
159 
160 /*
161  class exchangeInfo{
162  PS::ReallocatableArray<PS::S32> n_send_;
163  PS::ReallocatableArray<PS::S32> n_recv_;
164  PS::ReallocatableArray<PS::S32> n_disp_send_;
165  PS::ReallocatableArray<PS::S32> n_disp_recv_;
166  void set(const PS::S32 n_proc_send, const PS::S32 n_proc_recv){
167  n_send_.resizeNoInitialize(n_proc_send);
168  n_disp_send_.resizeNoInitialize(n_proc_send+1);
169  n_recv_.resizeNoInitialize(n_proc_recv);
170  n_disp_recv_.resizeNoInitialize(n_proc_recv+1);
171  }
172  };
173 */
174 
175 
176 class SearchCluster{
177 public:
178  PS::ReallocatableArray<PS::S32> n_ptcl_in_multi_cluster_isolated_;
179  PS::ReallocatableArray<PS::S32> n_ptcl_in_multi_cluster_isolated_offset_;
180  PS::ReallocatableArray<Mediator> mediator_sorted_id_cluster_;
181  PS::ReallocatableArray<PtclComm> ptcl_recv_;
182  PS::ReallocatableArray<PS::S32> adr_sys_multi_cluster_isolated_;
183 private:
184  PS::ReallocatableArray<Cluster> cluster_comm_;
185  std::unordered_map<PS::S32, PS::S32> id_to_adr_pcluster_;
186  // 1st: adr of self 2nd: adr of ngb
187  PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > adr_ngb_multi_cluster_;
188  PS::ReallocatableArray<PS::S32> * adr_sys_one_cluster_;
189  PS::ReallocatableArray<PtclCluster> * ptcl_cluster_;
190  PS::ReallocatableArray<PS::S32> id_cluster_send_;
191  PS::ReallocatableArray<PS::S32> id_cluster_recv_;
192  PS::ReallocatableArray<PS::S32> adr_pcluster_send_;
193  PS::ReallocatableArray<PS::S32> adr_pcluster_recv_;
194  PS::ReallocatableArray<PS::S32> n_cluster_send_;
195  PS::ReallocatableArray<PS::S32> n_cluster_recv_;
196  PS::ReallocatableArray<PS::S32> n_cluster_disp_send_;
197  PS::ReallocatableArray<PS::S32> n_cluster_disp_recv_;
198  PS::ReallocatableArray<PS::S32> rank_send_cluster_;
199  PS::ReallocatableArray<PS::S32> rank_recv_cluster_;
200  PS::S32 n_pcluster_self_node_;
201  PS::ReallocatableArray<PS::S32> adr_sys_ptcl_send_;
202  PS::ReallocatableArray<PtclComm> ptcl_send_;
203  PS::ReallocatableArray<PS::S32> rank_send_ptcl_;
204  PS::ReallocatableArray<PS::S32> n_ptcl_send_;
205  PS::ReallocatableArray<PS::S32> n_ptcl_disp_send_;
206  PS::ReallocatableArray<PS::S32> rank_recv_ptcl_;
207  PS::ReallocatableArray<PS::S32> n_ptcl_recv_;
208  PS::ReallocatableArray<PS::S32> n_ptcl_disp_recv_;
209  template<class T>
210  void packDataToThread0(T * data){
211  const PS::S32 n_thread = PS::Comm::getNumberOfThread();
212  PS::S32 size_data_th0 = data[0].size();
213  PS::S32 size_data = size_data_th0;
214  for(PS::S32 i=1; i<n_thread; i++){
215  size_data += data[i].size();
216  }
217 #ifdef HARD_DEBUG
218  assert(size_data<ARRAY_ALLOW_LIMIT);
219 #endif
220 
221  data[0].resizeNoInitialize(size_data);
222 #pragma omp parallel
223  {
224  const PS::S32 ith = PS::Comm::getThreadNum();
225  if(ith > 0){
226  PS::S32 offset = size_data_th0;
227  for(PS::S32 i=1; i<ith; i++) offset += data[i].size();
228  for(PS::S32 i=0; i<data[ith].size(); i++) data[0][i+offset] = data[ith][i];
229  }
230  }
231  }
232 
233  void searchClusterImpl(PtclCluster * target,
234  PtclCluster * p_first,
235  PtclCluster *& p_top,
236  PS::S32 & n_ptcl_in_cluster,
237  PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > & adr_ngb_array,
238  bool & fg_isolated){
239  if(target->adr_sys_ < 0) fg_isolated = false;
240  target->flag_searched_ = true;
241  target->next_ = NULL;
242  n_ptcl_in_cluster++;
243  PS::S32 target_adr = target->adr_ngb_head_;
244  const PS::S32 n_ngb = target->n_ngb_;
245  for(PS::S32 i=0; i<n_ngb; i++){
246  PS::S32 adr_ngb = adr_ngb_array[target_adr+i].second;
247  if( (p_first+adr_ngb)->flag_searched_ == false){
248  p_top->next_ = p_first+adr_ngb;
249  p_top = p_first+adr_ngb;
250  searchClusterImpl(p_first+adr_ngb, p_first, p_top, n_ptcl_in_cluster,
251  adr_ngb_array, fg_isolated);
252  }
253  }
254  }
255 
256  void setNgbAdrHead(PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > id_ngb_multi_cluster[]){
257  const PS::S32 size_ngb = id_ngb_multi_cluster[0].size();
258  PS::S32 n_cnt = 0;
259  for(PS::S32 i=0; i<size_ngb; i++){
260  if(id_ngb_multi_cluster[0][i].first == ptcl_cluster_[0][n_cnt].id_){
261  ptcl_cluster_[0][n_cnt].adr_ngb_head_ = i;
262  n_cnt++;
263  }
264  }
265  }
266 
267  void mergePtclCluster(const PS::ReallocatableArray<PtclOuter> ptcl_outer[],
268  PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > id_ngb_multi_cluster[]){
269  if(ptcl_outer[0].size() == 0) return;
270  const PS::S32 size = ptcl_outer[0].size();
271  PS::S32 id_tmp = ptcl_outer[0][0].id_;
272  PS::S32 adr_ngb_head = id_ngb_multi_cluster[0].size();
273  PS::S32 rank_tmp = ptcl_outer[0][0].rank_org_;
274  PS::S32 n_cnt = 0;
275  for(PS::S32 i=0; i<size; i++){
276  if(id_tmp != ptcl_outer[0][i].id_){
277  ptcl_cluster_[0].push_back( PtclCluster(id_tmp, -1, adr_ngb_head, n_cnt, false, NULL, rank_tmp) );
278  id_tmp = ptcl_outer[0][i].id_;
279  rank_tmp = ptcl_outer[0][i].rank_org_;
280  adr_ngb_head += n_cnt;
281  n_cnt = 0;
282  }
283  id_ngb_multi_cluster[0].push_back( std::pair<PS::S32, PS::S32>
284  (ptcl_outer[0][i].id_,
285  ptcl_outer[0][i].id_ngb_));
286  n_cnt++;
287  }
288  ptcl_cluster_[0].push_back( PtclCluster(id_tmp, -1, adr_ngb_head, n_cnt, false, NULL, rank_tmp) );
289  }
290 
291  struct OPEqualID{
292  template<class T> bool operator() (const T & left, const T & right) const {
293  return left.id_ == right.id_;
294  }
295  };
296  struct OPEqualSecond{
297  template<class T> bool operator() (const T & left, const T & right) const {
298  return left.second == right.second;
299  }
300  };
301  struct OPLessID{
302  template<class T> bool operator() (const T & left, const T & right) const {
303  return left.id_ < right.id_;
304  }
305  };
306  struct OPLessIDCluster{
307  template<class T> bool operator() (const T & left, const T & right) const {
308  return left.id_cluster_ < right.id_cluster_;
309  }
310  };
311  struct OPLessRankOrg{
312  template<class T> bool operator() (const T & left, const T & right) const {
313  return left.rank_org_ < right.rank_org_;
314  }
315  };
316  struct OPLessRank{
317  template<class T> bool operator() (const T & left, const T & right) const {
318  return left.rank_ < right.rank_;
319  }
320  };
321  struct OPLessFirst{
322  template<class T> bool operator() (const T & left, const T & right) const {
323  return left.first < right.first;
324  }
325  };
326  struct OPLessSecond{
327  template<class T> bool operator() (const T & left, const T & right) const {
328  return left.second < right.second;
329  }
330  };
331 
332 public:
333  void initialize(){
334  const PS::S32 n_thread = PS::Comm::getNumberOfThread();
335  adr_sys_one_cluster_ = new PS::ReallocatableArray<PS::S32>[n_thread];
336  ptcl_cluster_ = new PS::ReallocatableArray<PtclCluster>[n_thread];
337  }
338 
339 
341 
350  template<class Tpsoft, class Tepj>
351  PS::S32 checkNeighborWithVelocity(PS::S32* _index, Tpsoft& _pi, Tepj *_pb, const PS::S32 _nb, const PS::F64 _G, const PS::F64 _radius_factor) {
352 
353  PS::S32 n_nb_new = 0;
354  PS::F64 r_crit_i = _pi.changeover.getRout();
355 
356  ParticleBase pi;
357  pi.DataCopy(_pi);
358  auto& pi_cm = _pi.group_data.cm;
359 
360 #ifdef CLUSTER_DEBUG
361  assert(pi_cm.mass>=0.0);
362 #endif
363  // if i particle is member
364  if (pi_cm.mass>0.0) {
365  pi.mass = pi_cm.mass;
366  pi.vel[0] = pi_cm.vel.x;
367  pi.vel[1] = pi_cm.vel.y;
368  pi.vel[2] = pi_cm.vel.z;
369 //#ifdef CLUSTER_DEBUG
370 // if (_pi.group_data.artificial.status==0.0 && _pi.group_data.artificial.mass_backup!=0.0)
371 // std::cout<<"Warning: may not be pcm data! idi "<<_pi.id<<" status="<<_pi.group_data.artificial.status<<" mass_bk="<<_pi.group_data.artificial.mass_backup<<std::endl;
372 //#endif
373  }
374 
375  for (PS::S32 j=0; j<_nb; j++) {
376  if (_pi.id==_pb[j].id) continue;
377  auto& pbj_cm = _pb[j].group_data.cm;
378 #ifdef CLUSTER_DEBUG
379  assert(pbj_cm.mass>=0.0);
380 #endif
381  ParticleBase pj;
382  pj.pos = _pb[j].pos;
383  // particle member case
384  if (pbj_cm.mass>0.0) {
385  pj.vel[0] = pbj_cm.vel.x;
386  pj.vel[1] = pbj_cm.vel.y;
387  pj.vel[2] = pbj_cm.vel.z;
388  pj.mass = pbj_cm.mass;
389 //#ifdef CLUSTER_DEBUG
390 // if (_pb[j].group_data.artificial.status==0.0 && _pb[j].group_data.artificial.mass_backup!=0.0)
391 // std::cout<<"Warning: may not be pcm data! idj "<<_pb[j].id<<" status="<<_pb[j].group_data.artificial.status<<" mass_bk="<<_pb[j].group_data.artificial.mass_backup<<std::endl;
392 //#endif
393  }
394  else {
395  pj.mass= _pb[j].mass;
396  pj.vel = _pb[j].vel;
397  }
398 
399  PS::F64 semi,ecc, r,rv;
400  COMM::Binary::particleToSemiEcc(semi, ecc, r, rv, pi, pj, _G);
401  PS::F64 peri = semi*(1-ecc);
402  PS::F64 r_crit_j = _pb[j].r_out;
403  PS::F64 r_crit_max = _radius_factor*std::max(r_crit_i, r_crit_j);
404 
405  if (peri < r_crit_max && !(rv<0 && r>r_crit_max) ) {
406  _index[n_nb_new++] = j;
407  }
408 #ifdef CLUSTER_DEBUG_PRINT
409  else
410  std::cerr<<"Reject id.i "<<_pi.id<<" id.j "<<_pb[j].id<<" peri "
411  <<peri<<" r_out.i "<<r_crit_i<<" r_out.j "<<r_crit_j<<" dr "<<r<<" drdv "<<rv
412  <<" cm.vel.i "<<pi.vel[0]<<" "<<pi.vel[1]<<" "<<pi.vel[2]<<" m.i "<<pi.mass<<" "
413  <<" status.i "<<_pi.group_data.artificial.status<<" mass_bk.i "<<_pi.group_data.artificial.mass_backup
414  <<" cm.vel.j "<<pj.vel[0]<<" "<<pj.vel[1]<<" "<<pj.vel[2]<<" m.j "<<pj.mass<<" "
415  <<" status.j "<<_pb[j].group_data.artificial.status<<" mass_bk.j "<<_pb[j].group_data.artificial.mass_backup
416  <<std::endl;
417 #endif
418  }
419 
420  return n_nb_new;
421  }
422 
424  template<class Tsys, class Ttree, class Tepj>
425  void searchNeighborOMP(Tsys & sys,
426  Ttree & tree,
427  const PS::F64ort pos_domain[],
428  const PS::F64 _G,
429  const PS::F64 _radius_factor){
430  static PS::ReallocatableArray<PtclOuter> * ptcl_outer = NULL;
431  static PS::ReallocatableArray< std::pair<PS::S32, PS::S32> > * id_ngb_multi_cluster = NULL;
432  const PS::S32 n_thread = PS::Comm::getNumberOfThread();
433  if(ptcl_outer==NULL) ptcl_outer = new PS::ReallocatableArray<PtclOuter>[n_thread];
434  if(id_ngb_multi_cluster==NULL) id_ngb_multi_cluster = new PS::ReallocatableArray< std::pair<PS::S32, PS::S32> >[n_thread];
435  const PS::S32 my_rank = PS::Comm::getRank();
436  // const PS::S32 n_proc_tot = PS::Comm::getNumberOfProc();
437  const PS::S32 n_loc = sys.getNumberOfParticleLocal();
438 #ifdef CLUSTER_VELOCITY
439  assert(Ptcl::group_data_mode==GroupDataMode::cm);
440 #endif
441 
442 #pragma omp parallel
443  {
444  const PS::S32 ith = PS::Comm::getThreadNum();
445 // const PS::S32 ith = 0;
446  adr_sys_one_cluster_[ith].clearSize();
447  id_ngb_multi_cluster[ith].clearSize();
448  for(PS::S32 i=0; i<ptcl_cluster_[ith].size(); i++) ptcl_cluster_[ith][i].clear();
449  ptcl_cluster_[ith].clearSize();
450  for(PS::S32 i=0; i<ptcl_outer[ith].size(); i++) ptcl_outer[ith][i].clear();
451  ptcl_outer[ith].clearSize();
452 #pragma omp for
453  for(PS::S32 i=0; i<n_loc; i++){
454  if(sys[i].n_ngb == 1){
455  // no neighbor
456  adr_sys_one_cluster_[ith].push_back(i);
457 #ifdef CLUSTER_DEBUG
458 // assert(sys[i].group_data.artificial.isSingle());
459 
460  Tepj * nbl = NULL;
461  PS::S32 n_ngb_tree_i = tree.getNeighborListOneParticle(sys[i], nbl) ;
462  if(n_ngb_tree_i!=sys[i].n_ngb) {
463  std::cerr<<"Error: particle "<<i<<" Tree neighbor search number ("<<n_ngb_tree_i<<") is inconsistent with force kernel neighbor number ("<<sys[i].n_ngb<<")!\n Neighbor ID from tree: ";
464  for(int j=0; j<n_ngb_tree_i; j++) {
465  std::cerr<<(nbl+j)->id<<" ";
466  }
467  std::cerr<<std::endl;
468  abort();
469  }
470 #endif
471  }
472  else{
473  // has neighbor
474  Tepj * nbl = NULL;
475 
476 #ifdef SAVE_NEIGHBOR_ID_IN_FORCE_KERNEL
477  Tepj nb_pack[4];
478  // use saved index to get epj
479  PS::S32 n_ngb_force_i = sys[i].n_ngb;
480 
481  if (n_ngb_force_i<=4) {
482  for (PS::S32 k=0; k<n_ngb_force_i; k++) {
483  nb_pack[k] = *tree.getEpjFromId(sys[i].id_ngb[k]);
484  }
485  nbl = nb_pack;
486  sys[i].n_ngb--;
487  }
488  else sys[i].n_ngb = tree.getNeighborListOneParticle(sys[i], nbl) - 1;
489 #else
490 #ifdef CLUSTER_DEBUG
491  PS::S32 n_ngb_force_i = sys[i].n_ngb;
492 #endif
493  sys[i].n_ngb = tree.getNeighborListOneParticle(sys[i], nbl) - 1;
494 #endif
495 #ifdef CLUSTER_DEBUG
496  if(n_ngb_force_i<sys[i].n_ngb+1) {
497  std::cerr<<"Error: particle "<<i<<" Tree neighbor search number ("<<sys[i].n_ngb+1<<") is inconsistent with force kernel neighbor number ("<<n_ngb_force_i<<")!\n Neighbor ID from tree: ";
498  for(int j=0; j<sys[i].n_ngb+1; j++) {
499  std::cerr<<(nbl+j)->id<<" ";
500  }
501  std::cerr<<std::endl;
502  abort();
503  }
504 #endif
505 
506  // no neighbor
507  if(sys[i].n_ngb == 0){
508 //#ifdef CLUSTER_DEBUG
509 // assert(sys[i].group_data.artificial.isSingle()); // in cm mode, this should not be used
510 //#endif
511  adr_sys_one_cluster_[ith].push_back(i);
512  continue;
513  }
514 
515  PS::S32 adr_ngb_head_i = id_ngb_multi_cluster[ith].size();
516  PS::S32 n_ngb_i = sys[i].n_ngb;
517 
518 #ifdef CLUSTER_VELOCITY
519  // Use velocity criterion to select neighbors
520  PS::S32 neighbor_index[n_ngb_i];
521  PS::S32 n_ngb_check = checkNeighborWithVelocity(neighbor_index, sys[i], nbl, n_ngb_i+1, _G, _radius_factor);
522 #ifdef CLUSTER_DEBUG
523  assert(n_ngb_check>=0);
524 #endif
525  sys[i].n_ngb = n_ngb_check;
526 
527  // no neighbor
528  if (n_ngb_check==0) {
529  adr_sys_one_cluster_[ith].push_back(i);
530  continue;
531  }
532  else {
533  for(PS::S32 j=0; j<n_ngb_check; j++) {
534  PS::S32 index = neighbor_index[j];
535  auto* nbj = nbl + index;
536 
537  id_ngb_multi_cluster[ith].push_back( std::pair<PS::S32, PS::S32>(sys[i].id, nbj->id) );
538  if( nbj->rank_org != my_rank ){
539  ptcl_outer[ith].push_back(PtclOuter(nbj->id, sys[i].id, nbj->rank_org));
540  }
541  }
542  ptcl_cluster_[ith].push_back( PtclCluster(sys[i].id, i, adr_ngb_head_i, n_ngb_check, false, NULL, my_rank) );
543  }
544 #else
545  // use neighbor lists from tree neighbor search
546  for (int j=0; j<n_ngb_i+1; j++) {
547  auto* nbj = nbl + j;
548  if (sys[i].id==nbj->id) continue;
549 
550  id_ngb_multi_cluster[ith].push_back( std::pair<PS::S32, PS::S32>(sys[i].id, nbj->id) );
551  if( nbj->rank_org != my_rank ){
552  ptcl_outer[ith].push_back(PtclOuter(nbj->id, sys[i].id, nbj->rank_org));
553  }
554  }
555  ptcl_cluster_[ith].push_back( PtclCluster(sys[i].id, i, adr_ngb_head_i, n_ngb_i, false, NULL, my_rank) );
556 #endif
557  }
558  }
559  } // end of OMP parallel
560  packDataToThread0(adr_sys_one_cluster_);
561  packDataToThread0(ptcl_cluster_);
562  packDataToThread0(id_ngb_multi_cluster);
563  packDataToThread0(ptcl_outer);
564  setNgbAdrHead(id_ngb_multi_cluster);
565  n_pcluster_self_node_ = ptcl_cluster_[0].size();
566  std::sort(ptcl_outer[0].getPointer(), ptcl_outer[0].getPointer(ptcl_outer[0].size()), OPLessID());
567 
568 /* if(PS::Comm::getRank() == 0 && ptcl_outer[0].size() > 0){
569  PS::S32 n_tmp = 1;
570  PS::S32 n_tmp_tmp = 0;
571  //std::cerr<<"ptcl_outer[0].size()="<<ptcl_outer[0].size()<<std::endl;
572  PS::S32 ref_tmp = ptcl_outer[0][0].id_;
573  for(PS::S32 i=0; i<ptcl_outer[0].size(); i++){
574  if( ref_tmp != ptcl_outer[0][i].id_){
575  ref_tmp = ptcl_outer[0][i].id_;
576  n_tmp++;
577  n_tmp_tmp = 0;
578  }
579  n_tmp_tmp++;
580  }
581  }*/
582  mergePtclCluster(ptcl_outer, id_ngb_multi_cluster); // ptcl_cluster_ is complited
583  id_to_adr_pcluster_.clear(); // temporally add
584  for(PS::S32 i=0; i<ptcl_cluster_[0].size(); i++){
585  id_to_adr_pcluster_.insert(std::pair<PS::S32, PS::S32>(ptcl_cluster_[0][i].id_, i));
586  }
587  adr_ngb_multi_cluster_.clearSize();
588 #ifdef HARD_DEBUG
589  assert(id_ngb_multi_cluster[0].size()<ARRAY_ALLOW_LIMIT);
590 #endif
591  adr_ngb_multi_cluster_.resizeNoInitialize(id_ngb_multi_cluster[0].size());
592  for(PS::S32 i=0; i<id_ngb_multi_cluster[0].size(); i++){
593  const PS::S32 adr_self = id_to_adr_pcluster_[id_ngb_multi_cluster[0][i].first];
594  const PS::S32 adr_ngb = id_to_adr_pcluster_[id_ngb_multi_cluster[0][i].second];
595  adr_ngb_multi_cluster_[i] = std::pair<PS::S32, PS::S32>(adr_self, adr_ngb);
596  }
597  }
598 
599  template<class Tsys>
600  void checkPtclCluster(const Tsys & sys){
601  for(PS::S32 i=0; i<n_pcluster_self_node_; i++){
602  if(ptcl_cluster_[0][i].id_ != sys[ptcl_cluster_[0][i].adr_sys_].id){
603  std::cerr<<"i="<<i<<" ptcl_cluster_[0][i].adr_sys_="<<ptcl_cluster_[0][i].adr_sys_
604  <<" ptcl_cluster_[0][i].id_="<<ptcl_cluster_[0][i].id_
605  <<" sys[ptcl_cluster_[0][i].adr_sys_].id="<<sys[ptcl_cluster_[0][i].adr_sys_].id<<std::endl;
606  }
607  assert(ptcl_cluster_[0][i].id_ == sys[ptcl_cluster_[0][i].adr_sys_].id);
608  }
609  std::cerr<<"PASS: checkPtclCluster"<<std::endl;
610  }
611 
612 
613  void searchClusterLocal(){
614  const PS::S32 my_rank = PS::Comm::getRank();
615  const PS::S32 n_loc = ptcl_cluster_[0].size();
616 
617  adr_sys_multi_cluster_isolated_.clearSize();
618  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
619  mediator_sorted_id_cluster_[i].clear();
620  }
621  mediator_sorted_id_cluster_.clearSize();
622  static PS::ReallocatableArray<Cluster> cluster_isolated;
623  for(PS::S32 i=0; i<cluster_isolated.size(); i++) cluster_isolated[i].clear();
624  cluster_isolated.clearSize();
625  n_ptcl_in_multi_cluster_isolated_.clearSize();
626  n_ptcl_in_multi_cluster_isolated_offset_.clearSize();
627  n_ptcl_in_multi_cluster_isolated_offset_.push_back(0);
628 
629  for(PS::S32 i=0; i<n_loc; i++){
630  bool flag_isolated = true;
631  if(ptcl_cluster_[0][i].flag_searched_ == false){
632  PS::S32 n_ptcl_in_cluster = 0;
633  PtclCluster * p_top = ptcl_cluster_[0].getPointer(i);
634  searchClusterImpl(ptcl_cluster_[0].getPointer(i), ptcl_cluster_[0].getPointer(),
635  p_top, n_ptcl_in_cluster,
636  adr_ngb_multi_cluster_, flag_isolated);
637  if(flag_isolated){
638  PtclCluster * p_tmp = ptcl_cluster_[0].getPointer(i);
639  PS::S32 id_cluster = p_tmp->id_;
640  PS::S32 adr_sys_multi_cluster_isolated_head = adr_sys_multi_cluster_isolated_.size();
641  while(p_tmp != NULL){
642  if(p_tmp->id_ < id_cluster) id_cluster = p_tmp->id_;
643  adr_sys_multi_cluster_isolated_.push_back(p_tmp->adr_sys_);
644  p_tmp = p_tmp->next_;
645  }
646  cluster_isolated.push_back( Cluster(id_cluster, n_ptcl_in_cluster, n_ptcl_in_cluster, adr_sys_multi_cluster_isolated_head, my_rank) );
647  n_ptcl_in_multi_cluster_isolated_.push_back(n_ptcl_in_cluster);
648  n_ptcl_in_multi_cluster_isolated_offset_.push_back(n_ptcl_in_cluster+n_ptcl_in_multi_cluster_isolated_offset_.back());
649  }
650  else{
651  PtclCluster * p_tmp = ptcl_cluster_[0].getPointer(i);
652  PS::S32 id_cluster = p_tmp->id_;
653  while(p_tmp != NULL){
654  if(p_tmp->id_ < id_cluster) id_cluster = p_tmp->id_;
655  PS::S32 adr_pcluster = p_tmp-ptcl_cluster_[0].getPointer();
656  mediator_sorted_id_cluster_.push_back(Mediator(p_tmp->id_, p_tmp->adr_sys_, adr_pcluster, -1, PS::Comm::getRank()) ); // temporarily send_rank_ is my_rank
657  p_tmp = p_tmp->next_;
658  }
659  }
660  }
661  }
662  }
663 
664  template<class Tsys>
665  void checkMediator(const Tsys & sys){
666  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
667  assert(mediator_sorted_id_cluster_[i].id_ == ptcl_cluster_[0][mediator_sorted_id_cluster_[i].adr_pcluster_].id_);
668  }
669  std::cerr<<"PASS: checkMediator (1st)"<<std::endl;
670  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
671  if( mediator_sorted_id_cluster_[i].adr_sys_ == -1) continue;
672  assert(mediator_sorted_id_cluster_[i].id_ == sys[mediator_sorted_id_cluster_[i].adr_sys_].id);
673  }
674  std::cerr<<"PASS: checkMediator (2nd)"<<std::endl;
675  }
676 
677  void setIdClusterLocal(){
678  for(PS::S32 i=0; i<ptcl_cluster_[0].size(); i++){
679  ptcl_cluster_[0][i].setIdClusterImpl(ptcl_cluster_[0], adr_ngb_multi_cluster_);
680  }
681  }
682 
683 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
684  template <class Ttree>
685  void connectNodes(PS::F64ort pos_domain[], Ttree & tree){
686  const PS::S32 n_proc_tot = PS::Comm::getNumberOfProc();
687  const PS::S32 my_rank = PS::Comm::getRank();
688  static PS::ReallocatableArray<PS::S32> rank_neighbor;
689  static PS::ReallocatableArray<PS::S32> n_send_org;
690  static PS::ReallocatableArray<PS::S32> n_send_org_disp;
691  static PS::ReallocatableArray<PS::S32> id_send_org;
692  static PS::ReallocatableArray<PS::S32> n_send;
693  static PS::ReallocatableArray<PS::S32> n_recv;
694  static PS::ReallocatableArray<PS::S32> n_send_disp;
695  static PS::ReallocatableArray<PS::S32> n_recv_disp;
696  static PS::ReallocatableArray<PS::S32> id_send;
697  static PS::ReallocatableArray<PS::S32> id_recv;
698  rank_neighbor.clearSize();
699  n_send_org.resizeNoInitialize(n_proc_tot);
700  n_send_org_disp.resizeNoInitialize(n_proc_tot+1);
701  SetRankComm(pos_domain, tree, rank_neighbor);
702  for(PS::S32 i=0; i<n_proc_tot; i++) n_send_org[i] = 0;
703 
704  for(PS::S32 i=n_pcluster_self_node_; i<ptcl_cluster_[0].size(); i++){
705  assert(ptcl_cluster_[0][i].rank_org_ != my_rank);
706  n_send_org[ptcl_cluster_[0][i].rank_org_]++;
707  }
708 
709  n_send_org_disp[0] = 0;
710  static PS::ReallocatableArray<PS::S32> n_cnt_send;
711  n_cnt_send.resizeNoInitialize(n_proc_tot);
712  for(PS::S32 i=0; i<n_proc_tot; i++) {
713  n_send_org_disp[i+1] = n_send_org_disp[i] + n_send_org[i];
714  n_cnt_send[i] = 0;
715  }
716 
717  id_send_org.resizeNoInitialize(n_send_org_disp[n_proc_tot]);
718  for(PS::S32 i=n_pcluster_self_node_; i<ptcl_cluster_[0].size(); i++){
719  PS::S32 rank = ptcl_cluster_[0][i].rank_org_;
720  PS::S32 adr = n_send_org_disp[rank] + n_cnt_send[rank];
721  id_send_org[adr] = ptcl_cluster_[0][i].id_;
722  n_cnt_send[rank]++;
723  }
724 
725  const PS::S32 n_rank_send = rank_neighbor.size();
726  n_send.resizeNoInitialize(n_rank_send);
727  n_recv.resizeNoInitialize(n_rank_send);
728  n_send_disp.resizeNoInitialize(n_rank_send+1);
729  n_recv_disp.resizeNoInitialize(n_rank_send+1);
730  id_send.resizeNoInitialize(ptcl_cluster_[0].size());
731 
732  n_send_disp[0] = 0;
733  for (PS::S32 i=0; i<n_rank_send; i++) {
734  PS::S32 rank = rank_neighbor[i];
735  n_send[i] = n_send_org[rank];
736  if (n_send[i]!=n_cnt_send[rank]) {
737  std::cerr<<"n_send["<<i<<"]="<<n_send[i]<<" != n_cnt_send["<<rank<<"]="<<n_cnt_send[rank]<<std::endl;
738  abort();
739  }
740  n_send_disp[i+1] = n_send_disp[i] + n_send[i];
741  n_recv[i] = 0;
742 
743  for (PS::S32 j=0; j<n_send[i]; j++) {
744  id_send[n_send_disp[i]+j] = id_send_org[n_send_org_disp[rank]+j];
745  }
746  }
747 
748 #ifdef FDPS_COMM
749  PS::Comm::sendIrecv(n_send.getPointer(), rank_neighbor.getPointer(), 1, n_rank_send,
750  n_recv.getPointer(), rank_neighbor.getPointer(), 1, n_rank_send);
751 #else
752 
753  static PS::ReallocatableArray<MPI_Request> req_send;
754  static PS::ReallocatableArray<MPI_Request> req_recv;
755  static PS::ReallocatableArray<MPI_Status> stat_send;
756  static PS::ReallocatableArray<MPI_Status> stat_recv;
757  req_send.resizeNoInitialize(n_proc_tot);
758  req_recv.resizeNoInitialize(n_proc_tot);
759  stat_send.resizeNoInitialize(n_proc_tot);
760  stat_recv.resizeNoInitialize(n_proc_tot);
761 
762  for(PS::S32 i=0; i<rank_neighbor.size(); i++){
763  PS::S32 rank = rank_neighbor[i];
764  MPI_Isend(n_send.getPointer(i), 1, PS::GetDataType<PS::S32>(),
765  rank, 1837, MPI_COMM_WORLD, req_send.getPointer(i));
766  MPI_Irecv(n_recv.getPointer(i), 1, PS::GetDataType<PS::S32>(),
767  rank, 1837, MPI_COMM_WORLD, req_recv.getPointer(i));
768  }
769  MPI_Waitall(rank_neighbor.size(), req_send.getPointer(), stat_send.getPointer());
770  MPI_Waitall(rank_neighbor.size(), req_recv.getPointer(), stat_recv.getPointer());
771 #endif
772 
773  n_recv_disp[0] = 0;
774  for(PS::S32 i=0; i<n_rank_send; i++){
775  n_recv_disp[i+1] = n_recv_disp[i] + n_recv[i];
776  }
777  id_recv.resizeNoInitialize(n_recv_disp[n_rank_send]);
778 
779 #ifdef FDPS_COMM
780  PS::Comm::sendIrecvV(id_send.getPointer(), rank_neighbor.getPointer(), n_send.getPointer(), n_send_disp.getPointer(), n_rank_send,
781  id_recv.getPointer(), rank_neighbor.getPointer(), n_recv.getPointer(), n_recv_disp.getPointer(), n_rank_send);
782 #endif
783 
784  rank_send_cluster_.resizeNoInitialize(0);
785  rank_recv_cluster_.resizeNoInitialize(0);
786  n_cluster_send_.resizeNoInitialize(0);
787  n_cluster_recv_.resizeNoInitialize(0);
788  PS::S32 n_proc_send = 0;
789  PS::S32 n_proc_recv = 0;
790  for(PS::S32 i=0; i<n_rank_send; i++){
791  PS::S32 rank = rank_neighbor[i];
792  if(n_send[i] > 0){
793  rank_recv_cluster_.push_back(rank); // NOTE: recv not send
794  n_cluster_recv_.push_back(n_send[i]); // NOTE: recv not send
795 #ifndef FDPS_COMM
796  MPI_Isend(id_send.getPointer(n_send_disp[i]), n_send[i], PS::GetDataType<PS::S32>(),
797  rank, 1871, MPI_COMM_WORLD, req_send.getPointer(n_proc_recv));
798 #endif
799  n_proc_recv++;
800  }
801  if(n_recv[i] > 0){
802  rank_send_cluster_.push_back(rank); // NOTE: send
803  n_cluster_send_.push_back(n_recv[i]); // NOTE: send
804 #ifndef FDPS_COMM
805  MPI_Irecv(id_recv.getPointer(n_recv_disp[i]), n_recv[i], PS::GetDataType<PS::S32>(),
806  rank, 1871, MPI_COMM_WORLD, req_recv.getPointer(n_proc_send));
807 #endif
808  n_proc_send++;
809  }
810  }
811 #ifndef FDPS_COMM
812  MPI_Waitall(n_proc_recv, req_send.getPointer(), stat_send.getPointer());
813  MPI_Waitall(n_proc_send, req_recv.getPointer(), stat_recv.getPointer());
814 #endif
815  adr_pcluster_send_.clearSize();
816  adr_pcluster_recv_.clearSize();
817  for(PS::S32 i=0; i<n_recv_disp.back(); i++){
818 #ifdef CLUSTER_DEBUG
819  auto itr = id_to_adr_pcluster_.find(id_recv[i]);
820  assert(itr != id_to_adr_pcluster_.end());
821  assert(id_to_adr_pcluster_[id_recv[i]] < n_pcluster_self_node_);
822 #endif
823  adr_pcluster_send_.push_back(id_to_adr_pcluster_[id_recv[i]]);
824  }
825  for(PS::S32 i=0; i<n_send_disp.back(); i++){
826 #ifdef CLUSTER_DEBUG
827  auto itr = id_to_adr_pcluster_.find(id_send[i]);
828  assert(itr != id_to_adr_pcluster_.end());
829  /*
830  if(id_to_adr_pcluster_[id_send[i]] < n_pcluster_self_node_ && my_rank == 0){
831  std::cout<<"my_rank= "<<my_rank
832  <<" id_send[i]= "<<id_send[i]
833  <<" id_to_adr_pcluster_[id_send[i]]= "<<id_to_adr_pcluster_[id_send[i]]
834  <<" n_pcluster_self_node_= "<<n_pcluster_self_node_
835  <<std::endl;
836  }
837  */
838  assert(id_to_adr_pcluster_[id_send[i]] >= n_pcluster_self_node_);
839 #endif
840  adr_pcluster_recv_.push_back(id_to_adr_pcluster_[id_send[i]]);
841  }
842  n_cluster_disp_send_.resizeNoInitialize(n_proc_send+1);
843  n_cluster_disp_recv_.resizeNoInitialize(n_proc_recv+1);
844  n_cluster_disp_send_[0] = n_cluster_disp_recv_[0] = 0;
845  for(PS::S32 i=0; i<n_proc_send; i++) n_cluster_disp_send_[i+1] = n_cluster_disp_send_[i] + n_cluster_send_[i];
846  for(PS::S32 i=0; i<n_proc_recv; i++) n_cluster_disp_recv_[i+1] = n_cluster_disp_recv_[i] + n_cluster_recv_[i];
847  }
848 
849  void setIdClusterGlobalIteration(){
850  // const PS::S32 my_rank = PS::Comm::getRank();
851 #ifdef HARD_DEBUG
852  //assert(n_cluster_disp_send_.back()<ARRAY_ALLOW_LIMIT);
853  if(n_cluster_disp_send_.back()>ARRAY_ALLOW_LIMIT) {
854  std::cerr<<"Error: size overflow: rank: "<<PS::Comm::getRank()<<" n_cluster_disp_send_.back()="<<n_cluster_disp_send_.back()<<" size="<<n_cluster_disp_send_.size()<<std::endl;
855  }
856  if(n_cluster_disp_recv_.back()>ARRAY_ALLOW_LIMIT) {
857  std::cerr<<"Error: size overflow: rank: "<<PS::Comm::getRank()<<" n_cluster_disp_recv_.back()="<<n_cluster_disp_recv_.back()<<" size="<<n_cluster_disp_recv_.size()<<std::endl;
858  }
859 #endif
860  id_cluster_send_.resizeNoInitialize(n_cluster_disp_send_.back());
861  id_cluster_recv_.resizeNoInitialize(n_cluster_disp_recv_.back());
862 
863 #ifndef FDPS_COMM
864  const PS::S32 n_proc_tot = PS::Comm::getNumberOfProc();
865  static PS::ReallocatableArray<MPI_Request> req_send;
866  static PS::ReallocatableArray<MPI_Request> req_recv;
867  static PS::ReallocatableArray<MPI_Status> stat_send;
868  static PS::ReallocatableArray<MPI_Status> stat_recv;
869  req_send.resizeNoInitialize(n_proc_tot);
870  req_recv.resizeNoInitialize(n_proc_tot);
871  stat_send.resizeNoInitialize(n_proc_tot);
872  stat_recv.resizeNoInitialize(n_proc_tot);
873 #endif
874 
875  bool flag_itr_glb = true;
876  PS::S32 n_loop = 0;
877  while(flag_itr_glb){
878  flag_itr_glb = false;
879 
880  for(PS::S32 i=0; i<adr_pcluster_send_.size(); i++){
881  PS::S32 adr = adr_pcluster_send_[i];
882  assert(adr <= ptcl_cluster_[0].size());
883  id_cluster_send_[i] = ptcl_cluster_[0][adr].id_cluster_;
884  }
885 #ifdef FDPS_COMM
886  PS::Comm::sendIrecvV(id_cluster_send_.getPointer(), rank_send_cluster_.getPointer(), n_cluster_send_.getPointer(), n_cluster_disp_send_.getPointer(), rank_send_cluster_.size(),
887  id_cluster_recv_.getPointer(), rank_recv_cluster_.getPointer(), n_cluster_recv_.getPointer(), n_cluster_disp_recv_.getPointer(), rank_recv_cluster_.size());
888 #else
889  PS::S32 n_proc_send = 0;
890  PS::S32 n_proc_recv = 0;
891  for(PS::S32 i=0; i<rank_send_cluster_.size(); i++){
892  PS::S32 rank = rank_send_cluster_[i];
893  MPI_Isend(id_cluster_send_.getPointer(n_cluster_disp_send_[i]), n_cluster_send_[i], PS::GetDataType<PS::S32>(),
894  rank, 1947, MPI_COMM_WORLD, req_send.getPointer(n_proc_send));
895  n_proc_send++;
896  }
897  for(PS::S32 i=0; i<rank_recv_cluster_.size(); i++){
898  PS::S32 rank = rank_recv_cluster_[i];
899  MPI_Irecv(id_cluster_recv_.getPointer(n_cluster_disp_recv_[i]), n_cluster_recv_[i], PS::GetDataType<PS::S32>(),
900  rank, 1947, MPI_COMM_WORLD, req_recv.getPointer(n_proc_recv));
901  n_proc_recv++;
902  }
903  MPI_Waitall(n_proc_send, req_send.getPointer(), stat_send.getPointer());
904  MPI_Waitall(n_proc_recv, req_recv.getPointer(), stat_recv.getPointer());
905 #endif
906 
907  bool flag_itr_loc = false;
908  for(PS::S32 i=0; i<id_cluster_recv_.size(); i++){
909  PS::S32 adr = adr_pcluster_recv_[i];
910  if(id_cluster_recv_[i] < ptcl_cluster_[0][adr].id_cluster_){
911  ptcl_cluster_[0][adr].id_cluster_ = id_cluster_recv_[i];
912  flag_itr_loc |= true;
913  }
914  }
915  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
916  PS::S32 adr = mediator_sorted_id_cluster_[i].adr_pcluster_;
917  PS::S32 id_clucster_prev = ptcl_cluster_[0][adr].id_cluster_;
918  ptcl_cluster_[0][adr].setIdClusterImpl(ptcl_cluster_[0], adr_ngb_multi_cluster_);
919  if( id_clucster_prev != ptcl_cluster_[0][adr].id_cluster_){
920  flag_itr_loc |= true;
921  }
922  }
923  flag_itr_glb = PS::Comm::synchronizeConditionalBranchOR(flag_itr_loc);
924  n_loop++;
925  }
926  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
927  PS::S32 adr = mediator_sorted_id_cluster_[i].adr_pcluster_;
928  mediator_sorted_id_cluster_[i].id_cluster_ = ptcl_cluster_[0][adr].id_cluster_;
929  }
930  std::sort(mediator_sorted_id_cluster_.getPointer(0), mediator_sorted_id_cluster_.getPointer(mediator_sorted_id_cluster_.size()), OPLessIDCluster());
931  }
932 
933  template<class Tsys>
934  void sendAndRecvCluster(const Tsys & sys){
935  PS::S32 my_rank = PS::Comm::getRank();
936  PS::S32 n_proc = PS::Comm::getNumberOfProc();
937  static PS::ReallocatableArray<Cluster> cluster_loc;
938  cluster_loc.clearSize();
939  if(mediator_sorted_id_cluster_.size() > 0){
940  PS::S32 id_cluster_ref = mediator_sorted_id_cluster_[0].id_cluster_;
941  cluster_loc.push_back( Cluster(id_cluster_ref, 0, 0, 0, my_rank) );
942  for(PS::S32 i=0; i<mediator_sorted_id_cluster_.size(); i++){
943  if( id_cluster_ref != mediator_sorted_id_cluster_[i].id_cluster_){
944  id_cluster_ref = mediator_sorted_id_cluster_[i].id_cluster_;
945  cluster_loc.push_back( Cluster(id_cluster_ref, 0, 0, i, my_rank) );
946  }
947  if(mediator_sorted_id_cluster_[i].adr_sys_>=0) cluster_loc.back().n_ptcl_stored_++;
948  cluster_loc.back().n_ptcl_++;
949  }
950  }
951  //std::cerr<<"sendAndRecvCluster 0: "<<my_rank<<std::endl;
952 
953  static PS::ReallocatableArray<PS::S32> n_cluster_recv;
954  static PS::ReallocatableArray<PS::S32> n_cluster_recv_disp;
955  n_cluster_recv.resizeNoInitialize(n_proc);
956  n_cluster_recv_disp.resizeNoInitialize(n_proc+1);
957  PS::S32 n_cluster_tot_loc = cluster_loc.size();
958  PS::Comm::allGather(&n_cluster_tot_loc, 1, n_cluster_recv.getPointer());
959  n_cluster_recv_disp[0] = 0;
960  for(PS::S32 i=0; i<n_proc; i++){
961  n_cluster_recv_disp[i+1] = n_cluster_recv[i] + n_cluster_recv_disp[i];
962  }
963 
964  static PS::ReallocatableArray<Cluster> cluster_recv;
965  cluster_recv.resizeNoInitialize(n_cluster_recv_disp[n_proc]);
966  if(n_cluster_tot_loc > 0){
967  PS::Comm::allGatherV(cluster_loc.getPointer(), n_cluster_tot_loc,
968  cluster_recv.getPointer(),
969  n_cluster_recv.getPointer(), n_cluster_recv_disp.getPointer());
970  }
971  else{
972  Cluster tmp;
973  PS::Comm::allGatherV(&tmp, n_cluster_tot_loc,
974  cluster_recv.getPointer(),
975  n_cluster_recv.getPointer(), n_cluster_recv_disp.getPointer());
976  }
977 
978  //std::cerr<<"sendAndRecvCluster 1: "<<my_rank<<std::endl;
979 
980  // exchange cluster info
982 
983  for(PS::S32 i0=0; i0<n_cluster_tot_loc; i0++){
984  const PS::S32 id_cluster = cluster_loc[i0].id_;
985  const PS::S32 n_ptcl_in_cluster = cluster_loc[i0].n_ptcl_stored_;
986  PS::S32 rank_send_tmp = my_rank;
987  PS::S32 n_ptcl_max = n_ptcl_in_cluster;
988  for(PS::S32 i1=0; i1<n_proc; i1++){
989  if(i1 == my_rank) continue;
990  /*
991  if(my_rank == 0){
992  std::cout<<"i1="<<i1<<std::endl;
993  std::cout<<"n_cluster_recv_disp[i1]="<<n_cluster_recv_disp[i1]
994  <<" n_cluster_recv_disp[i1+1]="<<n_cluster_recv_disp[i1+1]<<std::endl;
995  }
996  */
997  for(PS::S32 i2=n_cluster_recv_disp[i1]; i2<n_cluster_recv_disp[i1+1]; i2++){
998  if(id_cluster == cluster_recv[i2].id_){
999  /*
1000  if(my_rank == 0){
1001  std::cout<<"cluster_recv[i2].id_="<<cluster_recv[i2].id_<<std::endl;
1002  std::cout<<"n_ptcl_in_cluster="<<n_ptcl_in_cluster<<std::endl;
1003  std::cout<<"cluster_recv[i2].n_ptcl_stored_="<<cluster_recv[i2].n_ptcl_stored_<<std::endl;
1004  }
1005  */
1006  if( (n_ptcl_max < cluster_recv[i2].n_ptcl_stored_) ||
1007  (n_ptcl_max == cluster_recv[i2].n_ptcl_stored_ && rank_send_tmp > i1) ){
1008  rank_send_tmp = i1;
1009  n_ptcl_max = cluster_recv[i2].n_ptcl_stored_;
1010  }
1011  }
1012  }
1013  }
1014  cluster_loc[i0].rank_ = rank_send_tmp;
1015  }
1016  std::sort(cluster_loc.getPointer(), cluster_loc.getPointer(cluster_loc.size()), OPLessRank());
1017  //std::cerr<<"sendAndRecvCluster 2: "<<my_rank<<std::endl;
1018 
1020  // pack and send particles
1021  ptcl_send_.clearSize();
1022  rank_send_ptcl_.clearSize();
1023  n_ptcl_send_.clearSize();
1024  adr_sys_ptcl_send_.clearSize();
1025  if(cluster_loc.size() > 0){
1026  PS::S32 rank_send_ref = -999999;
1027  for(PS::S32 i=0; i<cluster_loc.size(); i++){
1028  if(cluster_loc[i].rank_ == my_rank) continue;
1029  if( rank_send_ref != cluster_loc[i].rank_){
1030  rank_send_ref = cluster_loc[i].rank_;
1031  rank_send_ptcl_.push_back(rank_send_ref);
1032  n_ptcl_send_.push_back(0);
1033  }
1034 
1035  PS::S32 n_tmp = cluster_loc[i].n_ptcl_;
1036  PS::S32 n_cnt = 0;
1037  for(PS::S32 ii=0; ii<n_tmp; ii++){
1038  PS::S32 adr_sys = mediator_sorted_id_cluster_[cluster_loc[i].adr_head_+ii].adr_sys_;
1039  if(adr_sys < 0){
1040  continue;
1041  }
1042  mediator_sorted_id_cluster_[cluster_loc[i].adr_head_+ii].rank_send_ = rank_send_ref; // 2006.09.06
1043  const auto &p = sys[adr_sys];
1044  adr_sys_ptcl_send_.push_back(adr_sys);
1045  ptcl_send_.push_back(PtclComm(p));
1046  ptcl_send_.back().id_cluster = cluster_loc[i].id_;
1047 #ifdef HARD_DEBUG
1048  if(ptcl_send_.back().mass==0&&ptcl_send_.back().group_data.artificial.isUnused()) {
1049  std::cerr<<"Error! sending particle is unused! adr="<<adr_sys<<std::endl;
1050  abort();
1051  }
1052 #endif
1053  n_ptcl_send_.back()++;
1054  n_cnt++;
1055  }
1056  /*
1057  if(cluster_loc[i].n_ptcl_stored_ != n_cnt){
1058  std::cerr<<"cluster_loc[i].n_ptcl_stored_="<<cluster_loc[i].n_ptcl_stored_
1059  <<" n_cnt="<<n_cnt
1060  <<std::endl;
1061  }
1062  */
1063  assert(cluster_loc[i].n_ptcl_stored_ == n_cnt);
1064  }
1065  }
1066 
1067  //std::cerr<<"sendAndRecvCluster 3: "<<my_rank<<std::endl;
1068 
1069  //static PS::ReallocatableArray<PS::S32> n_ptcl_disp_send;
1070  n_ptcl_disp_send_.resizeNoInitialize(rank_send_ptcl_.size()+1);
1071  n_ptcl_disp_send_[0] = 0;
1072  for(PS::S32 i=0; i<rank_send_ptcl_.size(); i++){
1073  n_ptcl_disp_send_[i+1] = n_ptcl_disp_send_[i] + n_ptcl_send_[i];
1074  }
1075 
1076 #ifndef FDPS_COMM
1077  static PS::ReallocatableArray<MPI_Request> req_send;
1078  static PS::ReallocatableArray<MPI_Status> stat_send;
1079  req_send.resizeNoInitialize(n_proc);
1080  stat_send.resizeNoInitialize(n_proc);
1081  for(PS::S32 i=0; i<rank_send_ptcl_.size(); i++){
1082  PS::S32 rank = rank_send_ptcl_[i];
1083  MPI_Isend(ptcl_send_.getPointer(n_ptcl_disp_send_[i]), n_ptcl_send_[i],
1084  PS::GetDataType<PtclComm>(),
1085  rank, 2136, MPI_COMM_WORLD, req_send.getPointer(i));
1086  }
1087 #endif
1088 
1089  //std::cerr<<"sendAndRecvCluster 4: "<<my_rank<<std::endl;
1090  // pack and send particles
1092 
1094  // make and recv particles
1095  rank_recv_ptcl_.clearSize();
1096  n_ptcl_recv_.clearSize();
1097  for(PS::S32 i0=0; i0<n_proc; i0++){
1098  if(i0 == my_rank) continue;
1099  bool flag_recv = false;
1100  n_ptcl_recv_.push_back(0);
1101  rank_recv_ptcl_.push_back(i0);
1102  for(PS::S32 i1=n_cluster_recv_disp[i0]; i1<n_cluster_recv_disp[i0+1]; i1++){
1103  PS::S32 id_cluster = cluster_recv[i1].id_;
1104  for(PS::S32 i2=0; i2<cluster_loc.size(); i2++){
1105  if(id_cluster == cluster_loc[i2].id_ && cluster_loc[i2].rank_ == my_rank){
1106  n_ptcl_recv_.back() += cluster_recv[i1].n_ptcl_stored_;
1107  flag_recv = true;
1108  }
1109  }
1110  }
1111  if(!flag_recv){
1112  n_ptcl_recv_.resizeNoInitialize(n_ptcl_recv_.size()-1);
1113  rank_recv_ptcl_.resizeNoInitialize(rank_recv_ptcl_.size()-1);
1114  }
1115  }
1116 
1117  //static PS::ReallocatableArray<PS::S32> n_ptcl_disp_recv_;
1118  n_ptcl_disp_recv_.resizeNoInitialize(n_ptcl_recv_.size()+1);
1119  n_ptcl_disp_recv_[0] = 0;
1120  for(PS::S32 i=0; i<n_ptcl_recv_.size(); i++){
1121  n_ptcl_disp_recv_[i+1] = n_ptcl_disp_recv_[i] + n_ptcl_recv_[i];
1122  }
1123  //std::cerr<<"sendAndRecvCluster 4.5: "<<my_rank<<std::endl;
1124 
1125  ptcl_recv_.resizeNoInitialize(n_ptcl_disp_recv_[n_ptcl_recv_.size()]);
1126 
1127 #ifdef FDPS_COMM
1128  PS::Comm::sendIrecvV(ptcl_send_.getPointer(), rank_send_ptcl_.getPointer(), n_ptcl_send_.getPointer(), n_ptcl_disp_send_.getPointer(), rank_send_ptcl_.size(),
1129  ptcl_recv_.getPointer(), rank_recv_ptcl_.getPointer(), n_ptcl_recv_.getPointer(), n_ptcl_disp_recv_.getPointer(), rank_recv_ptcl_.size());
1130 
1131 #else
1132  static PS::ReallocatableArray<MPI_Request> req_recv;
1133  static PS::ReallocatableArray<MPI_Status> stat_recv;
1134  req_recv.resizeNoInitialize(n_proc);
1135  stat_recv.resizeNoInitialize(n_proc);
1136  for(PS::S32 i=0; i<rank_recv_ptcl_.size(); i++){
1137  PS::S32 rank = rank_recv_ptcl_[i];
1138  MPI_Irecv(ptcl_recv_.getPointer(n_ptcl_disp_recv_[i]), n_ptcl_recv_[i],
1139  PS::GetDataType<PtclComm>(),
1140  rank, 2136, MPI_COMM_WORLD, req_recv.getPointer(i));
1141  }
1142  //std::cerr<<"sendAndRecvCluster 4.7: "<<my_rank<<std::endl;
1143  MPI_Waitall(rank_send_ptcl_.size(), req_send.getPointer(), stat_send.getPointer());
1144  MPI_Waitall(rank_recv_ptcl_.size(), req_recv.getPointer(), stat_recv.getPointer());
1145 
1146  //std::cerr<<"sendAndRecvCluster 5: "<<my_rank<<std::endl;
1147 #endif
1148 
1149  // make and recv particles
1151  }
1152 
1154  /* Send local single particles to remote nodes and receive remote single particles.
1155  Notice _ptcl_hard are not overlap with ptcl_send
1156  @param[in,out] _sys: particle system. Notice the local particles of non-group members are updated.
1157  @param[in,out] _ptcl_hard: local partical in system_hard_connected
1158  */
1159  template<class Tsys, class Tphard>
1160  void SendSinglePtcl(Tsys & _sys,
1161  PS::ReallocatableArray<Tphard> & _ptcl_hard){
1162  for(PS::S32 i=0; i<ptcl_send_.size(); i++){
1163  PS::S32 adr = adr_sys_ptcl_send_[i];
1164 #ifdef HARD_DEBUG
1165  assert(_sys[adr].id == ptcl_send_[i].id);
1166 #endif
1167  // write local single particle
1168  if(ptcl_send_[i].group_data.artificial.isSingle() || (ptcl_send_[i].group_data.artificial.isMember() && ptcl_send_[i].getParticleCMAddress()<0)) {
1169  ptcl_send_[i].DataCopy(_sys[adr]);
1170  }
1171  }
1172 
1173 #ifdef FDPS_COMM
1174  PS::Comm::sendIrecvV(ptcl_send_.getPointer(), rank_send_ptcl_.getPointer(), n_ptcl_send_.getPointer(), n_ptcl_disp_send_.getPointer(), rank_send_ptcl_.size(),
1175  ptcl_recv_.getPointer(), rank_recv_ptcl_.getPointer(), n_ptcl_recv_.getPointer(), n_ptcl_disp_recv_.getPointer(), rank_recv_ptcl_.size());
1176 #else
1177 
1178  static PS::ReallocatableArray<MPI_Request> req_send;
1179  static PS::ReallocatableArray<MPI_Status> stat_send;
1180  static PS::ReallocatableArray<MPI_Request> req_recv;
1181  static PS::ReallocatableArray<MPI_Status> stat_recv;
1182  req_send.resizeNoInitialize(rank_send_ptcl_.size());
1183  stat_send.resizeNoInitialize(rank_send_ptcl_.size());
1184  req_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1185  stat_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1186 
1187  for(PS::S32 i=0; i<rank_send_ptcl_.size(); i++){
1188  PS::S32 rank = rank_send_ptcl_[i];
1189  MPI_Isend(ptcl_send_.getPointer(n_ptcl_disp_send_[i]), n_ptcl_send_[i],
1190  PS::GetDataType<PtclComm>(),
1191  rank, 2239, MPI_COMM_WORLD, req_send.getPointer(i));
1192  }
1193  for(PS::S32 i=0; i<rank_recv_ptcl_.size(); i++){
1194  PS::S32 rank = rank_recv_ptcl_[i];
1195  MPI_Irecv(ptcl_recv_.getPointer(n_ptcl_disp_recv_[i]), n_ptcl_recv_[i],
1196  PS::GetDataType<PtclComm>(),
1197  rank, 2239, MPI_COMM_WORLD, req_recv.getPointer(i));
1198  }
1199  MPI_Waitall(rank_send_ptcl_.size(), req_send.getPointer(), stat_send.getPointer());
1200  MPI_Waitall(rank_recv_ptcl_.size(), req_recv.getPointer(), stat_recv.getPointer());
1201 #endif
1202 
1203  // Receive remote single particle data
1204  const PS::S32 n = _ptcl_hard.size();
1205  for(PS::S32 i=0; i<n; i++){
1206  const PS::S32 adr = _ptcl_hard[i].adr_org;
1207  if(adr <0 && (_ptcl_hard[i].group_data.artificial.isSingle() || (_ptcl_hard[i].group_data.artificial.isMember() && _ptcl_hard[i].getParticleCMAddress()<0) )){
1208 #ifdef HARD_DEBUG
1209  assert( ptcl_recv_[-(adr+1)].id == _ptcl_hard[i].id );
1210 #endif
1211  _ptcl_hard[i].DataCopy(ptcl_recv_[-(adr+1)]);
1212  }
1213  }
1214 
1215  }
1216 
1218  /* Send local single particles to remote nodes and receive remote single particles.
1219  Notice _ptcl_hard are not overlap with ptcl_send
1220  @param[in,out] _sys: particle system.
1221  @param[in,out] _ptcl_hard: local partical in system_hard_connected
1222  */
1223  template<class Tsys, class Tphard>
1224  void SendPtcl(Tsys & _sys,
1225  PS::ReallocatableArray<Tphard> & _ptcl_hard){
1226  for(PS::S32 i=0; i<ptcl_send_.size(); i++){
1227  PS::S32 adr = adr_sys_ptcl_send_[i];
1228 #ifdef HARD_DEBUG
1229  assert(_sys[adr].id == ptcl_send_[i].id);
1230 #endif
1231  // write local single particle
1232  ptcl_send_[i].DataCopy(_sys[adr]);
1233  }
1234 
1235 #ifdef FDPS_COMM
1236  PS::Comm::sendIrecvV(ptcl_send_.getPointer(), rank_send_ptcl_.getPointer(), n_ptcl_send_.getPointer(), n_ptcl_disp_send_.getPointer(), rank_send_ptcl_.size(),
1237  ptcl_recv_.getPointer(), rank_recv_ptcl_.getPointer(), n_ptcl_recv_.getPointer(), n_ptcl_disp_recv_.getPointer(), rank_recv_ptcl_.size());
1238 #else
1239 
1240  static PS::ReallocatableArray<MPI_Request> req_send;
1241  static PS::ReallocatableArray<MPI_Status> stat_send;
1242  static PS::ReallocatableArray<MPI_Request> req_recv;
1243  static PS::ReallocatableArray<MPI_Status> stat_recv;
1244  req_send.resizeNoInitialize(rank_send_ptcl_.size());
1245  stat_send.resizeNoInitialize(rank_send_ptcl_.size());
1246  req_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1247  stat_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1248 
1249  for(PS::S32 i=0; i<rank_send_ptcl_.size(); i++){
1250  PS::S32 rank = rank_send_ptcl_[i];
1251  MPI_Isend(ptcl_send_.getPointer(n_ptcl_disp_send_[i]), n_ptcl_send_[i],
1252  PS::GetDataType<PtclComm>(),
1253  rank, 2239, MPI_COMM_WORLD, req_send.getPointer(i));
1254  }
1255  for(PS::S32 i=0; i<rank_recv_ptcl_.size(); i++){
1256  PS::S32 rank = rank_recv_ptcl_[i];
1257  MPI_Irecv(ptcl_recv_.getPointer(n_ptcl_disp_recv_[i]), n_ptcl_recv_[i],
1258  PS::GetDataType<PtclComm>(),
1259  rank, 2239, MPI_COMM_WORLD, req_recv.getPointer(i));
1260  }
1261  MPI_Waitall(rank_send_ptcl_.size(), req_send.getPointer(), stat_send.getPointer());
1262  MPI_Waitall(rank_recv_ptcl_.size(), req_recv.getPointer(), stat_recv.getPointer());
1263 #endif
1264 
1265  // Receive remote single particle data
1266  const PS::S32 n = _ptcl_hard.size();
1267  for(PS::S32 i=0; i<n; i++){
1268  const PS::S32 adr = _ptcl_hard[i].adr_org;
1269  if(adr <0){
1270 #ifdef HARD_DEBUG
1271  assert( ptcl_recv_[-(adr+1)].id == _ptcl_hard[i].id );
1272 #endif
1273  _ptcl_hard[i].DataCopy(ptcl_recv_[-(adr+1)]);
1274  }
1275  }
1276 
1277  }
1278 
1280  /* @param[in,out] _sys: particle system
1281  @param[in] _ptcl_hard: local partical in system_hard_connected
1282  @param[out] _mass_modify_list: address on _sys of particles where mass is modified
1283  */
1284  template<class Tsys, class Tphard>
1285  void writeAndSendBackPtcl(Tsys & _sys,
1286  const PS::ReallocatableArray<Tphard> & _ptcl_hard,
1287  PS::ReallocatableArray<PS::S32> &_mass_modify_list) {
1288  // const PS::S32 my_rank = PS::Comm::getRank();
1289  // const PS::S32 n_proc = PS::Comm::getNumberOfProc();
1290  const PS::S32 n = _ptcl_hard.size();
1291  for(PS::S32 i=0; i<n; i++){
1292  const PS::S32 adr = _ptcl_hard[i].adr_org;
1293  if( adr >= 0){
1294 #ifdef HARD_DEBUG
1295  assert( _sys[adr].id == _ptcl_hard[i].id);
1296 #endif
1297 #ifdef STELLAR_EVOLUTION
1298  PS::F64 mass_bk = _sys[adr].group_data.artificial.isMember()? _sys[adr].group_data.artificial.getMassBackup(): _sys[adr].mass;
1299  assert(mass_bk!=0.0);
1300 #endif
1301  _sys[adr].DataCopy(_ptcl_hard[i]);
1302 #ifdef STELLAR_EVOLUTION
1303  _sys[adr].dm = _sys[adr].mass - mass_bk;
1304  if (_sys[adr].dm!=0.0) _mass_modify_list.push_back(adr);
1305 #endif
1306  assert(!std::isinf(_sys[adr].pos[0]));
1307  assert(!std::isnan(_sys[adr].pos[0]));
1308  assert(!std::isinf(_sys[adr].vel[0]));
1309  assert(!std::isnan(_sys[adr].vel[0]));
1310  }
1311  else{
1312  //assert( ptcl_recv_[-(adr+1)].id == _ptcl_hard[i].id );
1313  ptcl_recv_[-(adr+1)].DataCopy(_ptcl_hard[i]);
1314  }
1315  }
1316 
1317 #ifdef FDPS_COMM
1318  PS::Comm::sendIrecvV(ptcl_recv_.getPointer(), rank_recv_ptcl_.getPointer(), n_ptcl_recv_.getPointer(), n_ptcl_disp_recv_.getPointer(), rank_recv_ptcl_.size(),
1319  ptcl_send_.getPointer(), rank_send_ptcl_.getPointer(), n_ptcl_send_.getPointer(), n_ptcl_disp_send_.getPointer(), rank_send_ptcl_.size());
1320 #else
1321 
1322  static PS::ReallocatableArray<MPI_Request> req_recv;
1323  static PS::ReallocatableArray<MPI_Status> stat_recv;
1324  static PS::ReallocatableArray<MPI_Request> req_send;
1325  static PS::ReallocatableArray<MPI_Status> stat_send;
1326  req_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1327  stat_recv.resizeNoInitialize(rank_recv_ptcl_.size());
1328  req_send.resizeNoInitialize(rank_send_ptcl_.size());
1329  stat_send.resizeNoInitialize(rank_send_ptcl_.size());
1330 
1331  for(PS::S32 i=0; i<rank_recv_ptcl_.size(); i++){
1332  PS::S32 rank = rank_recv_ptcl_[i];
1333  MPI_Isend(ptcl_recv_.getPointer(n_ptcl_disp_recv_[i]), n_ptcl_recv_[i],
1334  PS::GetDataType<PtclComm>(),
1335  rank, 2303, MPI_COMM_WORLD, req_recv.getPointer(i));
1336  }
1337  for(PS::S32 i=0; i<rank_send_ptcl_.size(); i++){
1338  PS::S32 rank = rank_send_ptcl_[i];
1339  MPI_Irecv(ptcl_send_.getPointer(n_ptcl_disp_send_[i]), n_ptcl_send_[i],
1340  PS::GetDataType<PtclComm>(),
1341  rank, 2303, MPI_COMM_WORLD, req_send.getPointer(i));
1342  }
1343  MPI_Waitall(rank_send_ptcl_.size(), req_send.getPointer(), stat_send.getPointer());
1344  MPI_Waitall(rank_recv_ptcl_.size(), req_recv.getPointer(), stat_recv.getPointer());
1345 #endif
1346 
1347  for(PS::S32 i=0; i<ptcl_send_.size(); i++){
1348  PS::S32 adr = adr_sys_ptcl_send_[i];
1349 #ifdef HARD_DEBUG
1350  assert(_sys[adr].id == ptcl_send_[i].id);
1351 #endif
1352 #ifdef STELLAR_EVOLUTION
1353  PS::F64 mass_bk = _sys[adr].group_data.artificial.isMember()? _sys[adr].group_data.artificial.getMassBackup(): _sys[adr].mass;
1354 #endif
1355  _sys[adr].DataCopy(ptcl_send_[i]);
1356 #ifdef STELLAR_EVOLUTION
1357  _sys[adr].dm = _sys[adr].mass - mass_bk;
1358  if (_sys[adr].dm!=0.0) _mass_modify_list.push_back(adr);
1359 #endif
1360  assert(!std::isinf(_sys[adr].pos[0]));
1361  assert(!std::isnan(_sys[adr].pos[0]));
1362  assert(!std::isinf(_sys[adr].vel[0]));
1363  assert(!std::isnan(_sys[adr].vel[0]));
1364  }
1365  }
1366 
1367 
1368 #endif
1369 
1370 
1371  const PS::ReallocatableArray<PS::S32> & getAdrSysOneCluster(){
1372  return adr_sys_one_cluster_[0];
1373  }
1374 
1376  const PS::ReallocatableArray<PS::S32> & getAdrSysConnectClusterSend(){
1377  return adr_sys_ptcl_send_;
1378  }
1379 
1380 };
1381 
1382 
Mediator::dump
void dump()
Definition: cluster_list.hpp:139
Cluster::n_ptcl_stored_
PS::S32 n_ptcl_stored_
Definition: cluster_list.hpp:53
PtclCluster::clear
void clear()
Definition: cluster_list.hpp:89
PtclCluster
Definition: cluster_list.hpp:72
PtclOuter
Definition: cluster_list.hpp:107
Cluster::n_ptcl_
PS::S32 n_ptcl_
Definition: cluster_list.hpp:52
PtclCluster::flag_searched_
bool flag_searched_
Definition: cluster_list.hpp:77
SearchCluster::checkPtclCluster
void checkPtclCluster(const Tsys &sys)
Definition: cluster_list.hpp:600
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
SearchCluster::n_ptcl_in_multi_cluster_isolated_
PS::ReallocatableArray< PS::S32 > n_ptcl_in_multi_cluster_isolated_
Definition: cluster_list.hpp:178
ParticleBase
Basic particle class.
Definition: particle_base.hpp:20
PtclCluster::next_
PtclCluster * next_
Definition: cluster_list.hpp:78
Ptcl::group_data_mode
static GroupDataMode group_data_mode
Definition: ptcl.hpp:47
PtclComm::PtclComm
PtclComm(const Tp &p)
Definition: cluster_list.hpp:151
SearchCluster::mediator_sorted_id_cluster_
PS::ReallocatableArray< Mediator > mediator_sorted_id_cluster_
Definition: cluster_list.hpp:180
PtclComm::print
void print(std::ostream &fout)
Definition: cluster_list.hpp:154
GroupDataMode::cm
@ cm
Mediator::rank_send_
PS::S32 rank_send_
Definition: cluster_list.hpp:129
Mediator::id_cluster_
PS::S32 id_cluster_
Definition: cluster_list.hpp:128
Cluster::Cluster
Cluster(const PS::S32 _id, const PS::S32 _n_ptcl, const PS::S32 _n_ptcl_stored, const PS::S32 _adr_head, const PS::S32 _rank)
Definition: cluster_list.hpp:57
SearchCluster::initialize
void initialize()
Definition: cluster_list.hpp:333
Mediator
Definition: cluster_list.hpp:123
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
Cluster
Definition: cluster_list.hpp:50
Cluster::adr_head_
PS::S32 adr_head_
Definition: cluster_list.hpp:54
PtclCluster::rank_org_
PS::S32 rank_org_
Definition: cluster_list.hpp:79
Cluster::rank_
PS::S32 rank_
Definition: cluster_list.hpp:55
ParticleBase::pos
PS::F64vec pos
Definition: particle_base.hpp:24
Ptcl
Particle class.
Definition: ptcl.hpp:36
ParticleBase::vel
PS::F64vec vel
Definition: particle_base.hpp:25
PtclOuter::PtclOuter
PtclOuter()
Definition: cluster_list.hpp:112
ptcl.hpp
PtclOuter::clear
void clear()
Definition: cluster_list.hpp:115
Cluster::id_
PS::S32 id_
Definition: cluster_list.hpp:51
SearchCluster::n_ptcl_in_multi_cluster_isolated_offset_
PS::ReallocatableArray< PS::S32 > n_ptcl_in_multi_cluster_isolated_offset_
Definition: cluster_list.hpp:179
PtclCluster::id_cluster_
PS::S32 id_cluster_
Definition: cluster_list.hpp:80
SearchCluster::checkMediator
void checkMediator(const Tsys &sys)
Definition: cluster_list.hpp:665
Mediator::clear
void clear()
Definition: cluster_list.hpp:136
PtclCluster::id_
PS::S32 id_
Definition: cluster_list.hpp:73
PtclOuter::id_ngb_
PS::S32 id_ngb_
Definition: cluster_list.hpp:110
Cluster::clear
void clear()
Definition: cluster_list.hpp:61
SearchCluster::getAdrSysConnectClusterSend
const PS::ReallocatableArray< PS::S32 > & getAdrSysConnectClusterSend()
get the address list of ptcl need to send in connected clusters
Definition: cluster_list.hpp:1376
PtclCluster::PtclCluster
PtclCluster()
Definition: cluster_list.hpp:81
PtclCluster::adr_ngb_head_
PS::S32 adr_ngb_head_
Definition: cluster_list.hpp:75
Cluster::dump
void dump()
Definition: cluster_list.hpp:64
SearchCluster::checkNeighborWithVelocity
PS::S32 checkNeighborWithVelocity(PS::S32 *_index, Tpsoft &_pi, Tepj *_pb, const PS::S32 _nb, const PS::F64 _G, const PS::F64 _radius_factor)
identify whether the neighbor satisfy velocity criterion
Definition: cluster_list.hpp:351
PtclCluster::n_ngb_
PS::S32 n_ngb_
Definition: cluster_list.hpp:76
SearchCluster::adr_sys_multi_cluster_isolated_
PS::ReallocatableArray< PS::S32 > adr_sys_multi_cluster_isolated_
Definition: cluster_list.hpp:182
Mediator::adr_pcluster_
PS::S32 adr_pcluster_
Definition: cluster_list.hpp:127
PtclOuter::id_
PS::S32 id_
Definition: cluster_list.hpp:109
SearchCluster::ptcl_recv_
PS::ReallocatableArray< PtclComm > ptcl_recv_
Definition: cluster_list.hpp:181
ParticleBase::DataCopy
void DataCopy(const Tp &din)
Copy from another ParticleBase.
Definition: particle_base.hpp:327
PtclOuter::dump
void dump()
Definition: cluster_list.hpp:118
PtclCluster::adr_sys_
PS::S32 adr_sys_
Definition: cluster_list.hpp:74
ARRAY_ALLOW_LIMIT
#define ARRAY_ALLOW_LIMIT
Definition: cluster_list.hpp:11
Mediator::Mediator
Mediator(const PS::S32 _id, const PS::S32 _adr_sys, const PS::S32 _adr_pcluster, const PS::S32 _id_cluster, const PS::S32 _rank_send)
Definition: cluster_list.hpp:131
SearchCluster::searchClusterLocal
void searchClusterLocal()
Definition: cluster_list.hpp:613
Mediator::Mediator
Mediator()
Definition: cluster_list.hpp:130
PtclOuter::rank_org_
PS::S32 rank_org_
Definition: cluster_list.hpp:111
ParticleBase::mass
PS::F64 mass
Definition: particle_base.hpp:23
PtclComm::id_cluster
PS::S32 id_cluster
Definition: cluster_list.hpp:148
PIKG::max
T max(const Vector3< T > &v)
Definition: pikg_vector.hpp:143
Cluster::Cluster
Cluster()
Definition: cluster_list.hpp:56
PtclComm
Definition: cluster_list.hpp:146
SearchCluster
Definition: cluster_list.hpp:176
Mediator::id_
PS::S32 id_
Definition: cluster_list.hpp:125
PtclOuter::PtclOuter
PtclOuter(const PS::S32 _id, const PS::S32 _id_ngb, const PS::S32 _rank_org)
Definition: cluster_list.hpp:113
PtclCluster::setIdClusterImpl
void setIdClusterImpl(PS::ReallocatableArray< PtclCluster > &ptcl, PS::ReallocatableArray< std::pair< PS::S32, PS::S32 > > &adr_ngb)
Definition: cluster_list.hpp:98
PtclComm::PtclComm
PtclComm()
Definition: cluster_list.hpp:152
Mediator::adr_sys_
PS::S32 adr_sys_
Definition: cluster_list.hpp:126
SearchCluster::setIdClusterLocal
void setIdClusterLocal()
Definition: cluster_list.hpp:677
SearchCluster::getAdrSysOneCluster
const PS::ReallocatableArray< PS::S32 > & getAdrSysOneCluster()
Definition: cluster_list.hpp:1371
SearchCluster::searchNeighborOMP
void searchNeighborOMP(Tsys &sys, Ttree &tree, const PS::F64ort pos_domain[], const PS::F64 _G, const PS::F64 _radius_factor)
search neighbors and separate isolated and multiple clusters
Definition: cluster_list.hpp:425
PtclCluster::dump
void dump()
Definition: cluster_list.hpp:93
PtclCluster::PtclCluster
PtclCluster(const PS::S32 _id, const PS::S32 _adr_sys, const PS::S32 _adr_ngb_head, const PS::S32 _n_ngb, const bool _flag_searched, PtclCluster *const _next, PS::S32 _rank_org)
Definition: cluster_list.hpp:83
Ptcl::print
void print(std::ostream &_fout) const
Definition: ptcl.hpp:72