PeTar
N-body code for collisional gravitational systems
domain.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 inline void CalculateBoundaryOfDomain(const PS::S32 &np,
4  const PS::F64vec pos_sample[],
5  const PS::S32 cid,
6  const PS::S32 &istart,
7  const PS::S32 &iend,
8  const PS::F64ort pos_root_domain,
9  PS::F64 & xlow,
10  PS::F64 & xhigh) {
11  if(istart == 0) {
12  xlow = pos_root_domain.low_[cid];
13  }
14  else {
15  xlow = 0.5 * (pos_sample[istart-1][cid] + pos_sample[istart][cid]);
16  }
17  if(iend == np - 1) {
18  xhigh = pos_root_domain.high_[cid];
19  }
20  else {
21  xhigh = 0.5 * (pos_sample[iend][cid] + pos_sample[iend+1][cid]);
22  }
23 }
24 
25 template<class Tsys>
26 inline void DomainDecision(PS::DomainInfo & dinfo,
27  const Tsys & system){
28  const PS::S32 n_proc = PS::Comm::getNumberOfProc();
29  assert(n_proc%4 == 0);
30  if(n_proc%4 == 0){
31  int n_loc_tmp = system.getNumberOfParticleLocal();
32  int * n_recv_tmp = new int[n_proc];
33  int * n_recv_disp_tmp = new int[n_proc+1];
34  PS::F64vec * pos_loc_tmp = new PS::F64vec[n_loc_tmp];
35  for(int i=0; i<n_loc_tmp; i++){
36  pos_loc_tmp[i].z = system[i].pos.z;
37  if(system[i].pos.x * system[i].pos.y > 0.0){
38  pos_loc_tmp[i].x = fabs(system[i].pos.x);
39  pos_loc_tmp[i].y = fabs(system[i].pos.y);
40  }
41  else{
42  pos_loc_tmp[i].x = fabs(system[i].pos.y);
43  pos_loc_tmp[i].y = fabs(system[i].pos.x);
44  }
45  }
46  PS::Comm::allGather(&n_loc_tmp, 1, n_recv_tmp);
47  n_recv_disp_tmp[0] = 0;
48  for(int i=0; i<n_proc; i++){
49  n_recv_disp_tmp[i+1] = n_recv_disp_tmp[i] + n_recv_tmp[i];
50  }
51  int n_glb_tmp = n_recv_disp_tmp[n_proc];
52  PS::F64vec * pos_glb_tmp = new PS::F64vec[n_glb_tmp];
53  PS::Comm::allGatherV(pos_loc_tmp, n_loc_tmp, pos_glb_tmp, n_recv_tmp, n_recv_disp_tmp);
54  PS::S32 n_proc_quat = n_proc / 4;
55  // fout_debug<<"n_proc, n_proc_quat="<<n_proc<<" "<<n_proc_quat<<std::endl;
56  PS::S32 nx_quat = sqrt((PS::F64)n_proc_quat-0.000001)+1;
57  while( n_proc_quat % nx_quat != 0) nx_quat++;
58  PS::S32 ny_quat = n_proc_quat / nx_quat;
59  PS::S32 nx = nx_quat*2;
60  PS::S32 ny = ny_quat*2;
61  PS::S32 nz = 1;
62 // if(PS::Comm::getRank() == 0){
63 // fout_debug<<"nx_quat, ny_quat, nx, ny, nz= "
64 // <<nx_quat<<" "<<ny_quat<<" "<<nx<<" "<<ny<<" "<<nz<<std::endl;
65 // }
66  dinfo.setNumberOfDomainMultiDimension(nx, ny, nz);
67  PS::F64ort * pos_domain_tmp = new PS::F64ort[n_proc];
68  if(PS::Comm::getRank() == 0){
69  // fout_debug<<"n_glb_tmp="<<n_glb_tmp<<std::endl;
70  PS::S32 * istart = new PS::S32[n_proc_quat];
71  PS::S32 * iend = new PS::S32[n_proc_quat];
72  PS::F64ort pos_root_domain_tmp = PS::F64ort( PS::F64vec(0.0, 0.0, -PS::LARGE_FLOAT), PS::F64vec(PS::LARGE_FLOAT, PS::LARGE_FLOAT, PS::LARGE_FLOAT));
73 #ifdef __HPC_ACE__
74  std::sort(pos_glb_tmp, pos_glb_tmp+n_glb_tmp, PS::LessOPX());
75 #else
76  std::sort(pos_glb_tmp, pos_glb_tmp+n_glb_tmp,
77  [](const PS::F64vec & l, const PS::F64vec & r)->bool{return l.x < r.x;}
78  );
79 #endif
80  // fout_debug<<"pos_glb_tmp[n_glb_tmp-1].x="<<pos_glb_tmp[n_glb_tmp-1].x<<std::endl;
81  for(PS::S32 i = 0; i < n_proc_quat; i++) {
82  istart[i] = ((PS::S64)(i) * (PS::S64)(n_glb_tmp)) / (PS::S64)(n_proc_quat);
83  if(i > 0) iend[i-1] = istart[i] - 1;
84  }
85  iend[n_proc_quat-1] = n_glb_tmp - 1;
86  for(PS::S32 ix = 0; ix<nx_quat; ix++) {
87  PS::S32 ix0 = ix * ny_quat * nz;
88  PS::S32 ix1 = (ix + 1) * ny_quat * nz;
89  PS::F64 x0, x1;
90  CalculateBoundaryOfDomain(n_glb_tmp, pos_glb_tmp, 0, istart[ix0], iend[ix1-1], pos_root_domain_tmp, x0, x1);
91  for(PS::S32 i=0; i<ny_quat*nz; i++) {
92  PS::S32 offset = (nx_quat+ix)*ny;
93  pos_domain_tmp[offset+i].low_[0] = x0;
94  pos_domain_tmp[offset+i].high_[0] = x1;
95  pos_domain_tmp[offset+ny_quat+i].low_[0] = x0;
96  pos_domain_tmp[offset+ny_quat+i].high_[0] = x1;
97 
98  offset = (nx_quat-1-ix)*ny;
99  pos_domain_tmp[offset+i].low_[0] = -x1;
100  pos_domain_tmp[offset+i].high_[0] = -x0;
101  pos_domain_tmp[offset+ny_quat+i].low_[0] = -x1;
102  pos_domain_tmp[offset+ny_quat+i].high_[0] = -x0;
103  }
104  }
105 
106  for(PS::S32 ix = 0; ix<nx_quat; ix++) {
107  PS::S32 ix0 = ix * ny_quat * nz;
108  PS::S32 ix1 = (ix + 1) * ny_quat * nz;
109 #ifdef __HPC_ACE__
110  std::sort(pos_glb_tmp+istart[ix0], pos_glb_tmp+iend[ix1-1]+1, PS::LessOPY());
111 #else
112  std::sort(pos_glb_tmp+istart[ix0], pos_glb_tmp+iend[ix1-1]+1,
113  [](const PS::F64vec & l, const PS::F64vec & r)->bool{return l.y < r.y;}
114  );
115 #endif
116  PS::S32 n_tmp_y = iend[ix1-1] - istart[ix0] + 1;
117  for(PS::S32 iy = 0; iy<ny_quat; iy++) {
118  PS::S32 iy0 = ix0 + iy * nz;
119  PS::S32 iy1 = ix0 + (iy + 1) * nz;
120  PS::F64 y0, y1;
121  CalculateBoundaryOfDomain(n_tmp_y, pos_glb_tmp+istart[ix0], 1, istart[iy0]-istart[ix0], iend[iy1-1]-istart[ix0], pos_root_domain_tmp, y0, y1);
122  for(PS::S32 i=0; i<nz; i++) {
123  PS::S32 offset = (nx_quat+ix)*ny + ny_quat + iy;
124  pos_domain_tmp[offset+i].low_[1] = y0;
125  pos_domain_tmp[offset+i].high_[1] = y1;
126  offset = (nx_quat-ix-1)*ny + ny_quat + iy;
127  pos_domain_tmp[offset+i].low_[1] = y0;
128  pos_domain_tmp[offset+i].high_[1] = y1;
129  offset = (nx_quat+ix)*ny + ny_quat-iy-1;
130  pos_domain_tmp[offset+i].low_[1] = -y1;
131  pos_domain_tmp[offset+i].high_[1] = -y0;
132  offset = (nx_quat-ix-1)*ny + ny_quat-iy-1;
133  pos_domain_tmp[offset+i].low_[1] = -y1;
134  pos_domain_tmp[offset+i].high_[1] = -y0;
135  }
136  }
137  }
138  delete [] istart;
139  delete [] iend;
140  }
141  PS::Comm::broadcast(pos_domain_tmp, n_proc);
142  for(PS::S32 i=0; i<n_proc; i++){
143  pos_domain_tmp[i].low_.z = -PS::LARGE_FLOAT;
144  pos_domain_tmp[i].high_.z = PS::LARGE_FLOAT;
145  dinfo.setPosDomain(i, pos_domain_tmp[i]);
146  }
147  delete [] n_recv_tmp;
148  delete [] n_recv_disp_tmp;
149  delete [] pos_loc_tmp;
150  delete [] pos_glb_tmp;
151  delete [] pos_domain_tmp;
152  }
153 }
DomainDecision
void DomainDecision(PS::DomainInfo &dinfo, const Tsys &system)
Definition: domain.hpp:26
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
PIKG::F64vec
Vector3< F64 > F64vec
Definition: pikg_vector.hpp:167
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
galpy_pot_movie.ny
ny
Definition: galpy_pot_movie.py:134
PIKG::S64
int64_t S64
Definition: pikg_vector.hpp:23
galpy_pot_movie.nz
nz
Definition: galpy_pot_movie.py:134
galpy_pot_movie.nx
nx
Definition: galpy_pot_movie.py:134
CalculateBoundaryOfDomain
void CalculateBoundaryOfDomain(const PS::S32 &np, const PS::F64vec pos_sample[], const PS::S32 cid, const PS::S32 &istart, const PS::S32 &iend, const PS::F64ort pos_root_domain, PS::F64 &xlow, PS::F64 &xhigh)
Definition: domain.hpp:3