PeTar
N-body code for collisional gravitational systems
search_group_candidate.hpp
Go to the documentation of this file.
1 #pragma once
2 #include<particle_simulator.hpp>
3 
4 #ifndef ARRAY_ALLOW_LIMIT
5 #define ARRAY_ALLOW_LIMIT 1000000000
6 #endif
7 
8 template<class Tptcl>
9 class SearchGroupCandidate{
10 private:
11  typedef std::pair<PS::S32, PS::S32> PLinker;
12 
13  // group member list
14  PS::ReallocatableArray<PS::S32> group_list_;
15  PS::ReallocatableArray<PS::S32> group_list_disp_;
16  PS::ReallocatableArray<PS::S32> group_list_n_;
17 
19  /* If r.v >0, Use bin factor sqrt(r^2 - (r*v/|v|)^2) to check, otherwise use distance to check. The mass ratio is also considered.
20  The perturber acceleration from non-partner member is stored in mass_bk for the later stablility check
21  @param[out] _part_list: partner index in _ptcl for each particle
22  @param[out] _part_list_disp: offset to separate partner index for different particles
23  @param[out] _part_list_n: number of partners for each particle
24  @param[in,out] _ptcl: particle data, mass_bk is updated
25  @param[in] _n: number of particles
26  */
27  void searchPartner(PS::ReallocatableArray<PS::S32> & _part_list,
28  PS::ReallocatableArray<PS::S32> & _part_list_disp,
29  PS::ReallocatableArray<PS::S32> & _part_list_n,
30  Tptcl *_ptcl,
31  const PS::S32 _n) {
32 
33  _part_list.clearSize();
34  _part_list_disp.reserve(_n);
35 #ifdef HARD_DEBUG
36  assert(_n<ARRAY_ALLOW_LIMIT);
37 #endif
38  _part_list_disp.resizeNoInitialize(_n);
39  _part_list_n.reserve(_n);
40  _part_list_n.resizeNoInitialize(_n);
41 
42  // find partner
43  PS::S32 offset = 0;
44  for(int i=0; i<_n; i++) {
45  //PS::S32 ip = p_list[i];
46  //_ptcl[i].mass_bk.d = 0.0;
47  _part_list_n[i] = 0;
48  _part_list_disp[i] = offset;
49  for(int j=0; j<_n; j++) {
50  //PS::S32 jp = p_list[j];
51  if(i==j) continue;
52  PS::F64vec dr = _ptcl[i].pos-_ptcl[j].pos;
53  PS::F64 r2 = dr*dr;
54  // use simple criterion
55  PS::F64 rin_min = std::min(_ptcl[i].getRGroupCandidate(), _ptcl[j].getRGroupCandidate());
56  if (r2<rin_min*rin_min) {
57  _part_list.push_back(j);
58  _part_list_n[i]++;
59  offset++;
60  }
61  else {
62  //store perturbation force
63  //_ptcl[i].mass_bk.d += _ptcl[j].mass/r2;
64  }
65 
66  //PS::F64vec dv = _ptcl[i].vel-_ptcl[j].vel;
67  //PS::F64 rv = dr*dv;
68  //PS::F64 r2 = dr*dr;
69  //PS::F64 mi = _ptcl[i].mass;
70  //PS::F64 mj = _ptcl[j].mass;
71  //PS::F64 mass_factor= mi>mj ? mi/mj : mj/mi; ///> get mass ratio
72  //PS::F64 mass_factor=std::max(ptcl[ip].mass,ptcl[jp].mass)*Ptcl::mean_mass_inv;
73 #ifdef HARD_DEBUG
74  if(r2==0) {
75  std::cerr<<"Error: zero distance! i="<<i<<" j="<<j<<std::endl;
76  abort();
77  }
78 #endif
79  // incoming case
80  //if (rv>0) {
81  // PS::F64 v2 = dv*dv;
82  // // use bin factor to check
83  // if (r2*v2-rv*rv<_r_crit2*v2*mass_factor) {
84  // _part_list.push_back(j);
85  // _part_list_n[i]++;
86  // offset++;
87  // }
88  // else {
89  // //store perturbation force
90  // _ptcl[i].mass_bk += _ptcl[j].mass/r2;
91  // }
92  //}
93  //else { // outgoing case
94  // if (r2<_r_crit2*mass_factor) {
95  // _part_list.push_back(j);
96  // _part_list_n[i]++;
97  // offset++;
98  // }
99  // else {
100  // //store perturbation force
101  // _ptcl[i].mass_bk += _ptcl[j].mass/r2;
102  // }
103  //}
104  }
105  }
106  }
107 
108  void mergeGroup(PS::ReallocatableArray<PS::S32> & group_list,
109  PS::ReallocatableArray<PS::S32> & group_list_disp,
110  PS::ReallocatableArray<PS::S32> & group_list_n,
111  //PS::ReallocatableArray<PS::S32> & p_list,
112  const PS::S32 _n_ptcl,
113  PS::S32 part_list[],
114  PS::S32 part_list_disp[],
115  PS::S32 part_list_n[]) {
116 
117  //const PS::S32 _n_ptcl = p_list.size();
118  // partner index with marker
119  PS::ReallocatableArray<PLinker> partner_index;
120  partner_index.reserve(_n_ptcl);
121 
122  // map index from ptcl_org to partner_index
123  PS::ReallocatableArray<PS::S32> reverse_list;
124 #ifdef HARD_DEBUG
125  assert(_n_ptcl<ARRAY_ALLOW_LIMIT);
126 #endif
127  reverse_list.reserve(_n_ptcl);
128  reverse_list.resizeNoInitialize(_n_ptcl);
129 
130 #ifdef HARD_DEBUG
131  assert(group_list.size()==0);
132  assert(group_list_disp.size()==0);
133  assert(group_list_n.size()==0);
134 #endif
135 
136  for(int i=0; i<_n_ptcl; i++) {
137  if(part_list_n[i]>0) {
138  reverse_list[i] = partner_index.size();
139  partner_index.push_back(PLinker(i,-1));
140  }
141 #ifdef HARD_DEBUG
142  else {
143  reverse_list[i] = -1;
144  }
145 #endif
146  }
147  PS::S32 n_tot = partner_index.size();
148  //PS::S32 n_groups = 0;
149 
150  PS::S32 n_mem = 0;
151  for(int i=0; i<n_tot; i++) {
152  // PS::S32 k = partner_index[i].first;
153  if(partner_index[i].second>=0) continue;
154 
155  PS::S32 npart = connectGroups(i,i,part_list, part_list_disp, part_list_n,partner_index,reverse_list);
156  group_list_n.push_back(npart);
157 #ifdef HARD_DEBUG
158  assert(npart>0);
159 #endif
160  group_list_disp_.push_back(n_mem);
161  n_mem += npart;
162  group_list.push_back(partner_index[i].first);
163  PS::S32 inext=partner_index[i].second;
164  PS::S32 k=1;
165  while (inext!=i) {
166  group_list.push_back(partner_index[inext].first);
167  inext=partner_index[inext].second;
168  k++;
169 #ifdef HARD_DEBUG
170  assert(k<=npart);
171 #endif
172  }
173 #ifdef HARD_DEBUG
174  assert(k==npart);
175  if(npart>_n_ptcl) {
176  std::cerr<<"Error: connect group particle number mismatch: npart ="<<npart<<" ; _n_ptcl = "<<_n_ptcl<<std::endl;
177  abort();
178  }
179 #endif
180 
181  //n_group.push_back(npart);
182  //n_groups++;
183  }
184  }
185 
186  PS::S32 connectGroups(const PS::S32 ip,
187  const PS::S32 iend,
188  PS::S32 part_list[],
189  PS::S32 part_list_disp[],
190  PS::S32 part_list_n[],
191  PS::ReallocatableArray<PLinker> & partner_index,
192  PS::ReallocatableArray<PS::S32> & reverse_list) {
193  PS::S32 n_connected = 0;
194  PS::S32 n_reduce = 0;
195  PS::U32 inow=ip;
196  PS::S32 kp = partner_index[ip].first;
197  std::vector<PS::U32> rlist;
198  for(PS::S32 j=0; j<part_list_n[kp]; j++) {
199  PS::S32 inext = reverse_list[part_list[part_list_disp[kp]+j]];
200  if(partner_index[inext].second<0&&inext!=iend) {
201  if(partner_index[inow].second>=0) n_reduce++;
202  partner_index[inow].second = inext;
203  rlist.push_back(inext);
204  n_connected++;
205  inow = inext;
206  }
207  }
208  if(n_connected>0) {
209  partner_index[inow].second = iend;
210  n_connected++;
211  }
212 
213  for(PS::U32 j=0; j<rlist.size(); j++) {
214  inow = rlist[j];
215  PS::U32 inext = partner_index[inow].second;
216  n_connected += connectGroups(inow,inext,part_list,part_list_disp,part_list_n,partner_index,reverse_list);
217  }
218  return n_connected - n_reduce;
219  }
220 
221 
222 public:
223 
224  void searchAndMerge(Tptcl *_ptcl_in_cluster, const PS::S32 _n_ptcl) {
225  PS::ReallocatableArray<PS::S32> part_list;
226  PS::ReallocatableArray<PS::S32> part_list_disp;
227  PS::ReallocatableArray<PS::S32> part_list_n;
228 
229  searchPartner(part_list, part_list_disp, part_list_n, _ptcl_in_cluster, _n_ptcl);
230  mergeGroup(group_list_, group_list_disp_, group_list_n_, _n_ptcl, part_list.getPointer(), part_list_disp.getPointer(), part_list_n.getPointer());
231  }
232 
233  PS::S32 getNumberOfGroups() const {
234  return group_list_n_.size();
235  }
236 
237  PS::S32* getMemberList(const std::size_t igroup) {
238  return &group_list_[group_list_disp_[igroup]];
239  }
240 
241  PS::S32 getGroupListSize() const {
242  return group_list_.size();
243  }
244 
245  PS::S32 getNumberOfGroupMembers(const std::size_t igroup) const {
246  return group_list_n_[igroup];
247  }
248 
249 };
250 
PIKG::U32
uint32_t U32
Definition: pikg_vector.hpp:21
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
PIKG::F64vec
Vector3< F64 > F64vec
Definition: pikg_vector.hpp:167
ARRAY_ALLOW_LIMIT
#define ARRAY_ALLOW_LIMIT
Definition: search_group_candidate.hpp:5
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
SearchGroupCandidate::getGroupListSize
PS::S32 getGroupListSize() const
Definition: search_group_candidate.hpp:241
SearchGroupCandidate::getMemberList
PS::S32 * getMemberList(const std::size_t igroup)
Definition: search_group_candidate.hpp:237
SearchGroupCandidate::getNumberOfGroups
PS::S32 getNumberOfGroups() const
Definition: search_group_candidate.hpp:233
SearchGroupCandidate::getNumberOfGroupMembers
PS::S32 getNumberOfGroupMembers(const std::size_t igroup) const
Definition: search_group_candidate.hpp:245
SearchGroupCandidate
Definition: search_group_candidate.hpp:9
SearchGroupCandidate::searchAndMerge
void searchAndMerge(Tptcl *_ptcl_in_cluster, const PS::S32 _n_ptcl)
Definition: search_group_candidate.hpp:224