PeTar
N-body code for collisional gravitational systems
status.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 #include "particle_base.hpp"
4 #include "energy.hpp"
5 
7 class Status {
8 public:
9  PS::F64 time; // time of system
10  PS::S64 n_real_loc; // number of real particles in the local MPI process
11  PS::S64 n_real_glb; // number of real particles in all MPI processes
12  PS::S64 n_all_loc; // number of all (real+aritficial) particles in the local MPI process
13  PS::S64 n_all_glb; // number of all (real+aritficial) particles in all MPI process
14  PS::S32 n_remove_glb; // number of removed particles in all MPI processes
15  PS::S32 n_escape_glb; // number of escaped particles in all MPI processes
16  PS::F64 half_mass_radius; // half mass radius of the system (not calculated)
17  EnergyAndMomemtum energy; // energy of the system
18  struct ParticleCM{ // local structure for system center
19  PS::F64 mass;
20  PS::F64vec pos;
21  PS::F64vec vel;
22  bool is_center_shift_flag;
23 
24  ParticleCM(): mass(0.0), pos(PS::F64vec(0.0)), vel(PS::F64vec(0.0)), is_center_shift_flag(false) {}
25 
27 
31  void printColumnTitle(std::ofstream & _fout, const PS::S32 _width=20) {
32  _fout<<std::setw(_width)<<"CM.mass"
33  <<std::setw(_width)<<"CM.pos.x"
34  <<std::setw(_width)<<"CM.pos.y"
35  <<std::setw(_width)<<"CM.pos.z"
36  <<std::setw(_width)<<"CM.vel.x"
37  <<std::setw(_width)<<"CM.vel.y"
38  <<std::setw(_width)<<"CM.vel.z";
39  }
40 
42 
46  void printColumn(std::ofstream & _fout, const PS::S32 _width=20) {
47  _fout<<std::setw(_width)<<mass
48  <<std::setw(_width)<<pos.x
49  <<std::setw(_width)<<pos.y
50  <<std::setw(_width)<<pos.z
51  <<std::setw(_width)<<vel.x
52  <<std::setw(_width)<<vel.y
53  <<std::setw(_width)<<vel.z;
54  }
55 
57 
60  void print(std::ostream & _fout) {
61  _fout<<"C.M.: mass: "<<mass
62  <<" pos: "<<pos
63  <<" vel: "<<vel;
64  }
65 
66  void clear() {
67  is_center_shift_flag = false;
68  mass = 0.0;
69  pos = PS::F64vec(0.0);
70  vel = PS::F64vec(0.0);
71  }
72 
73  } pcm;
74 
75  Status(): time(0.0), n_real_loc(0), n_real_glb(0), n_all_loc(0), n_all_glb(0), n_remove_glb(0), n_escape_glb(0), half_mass_radius(0), energy(), pcm() {}
76 
78 
82  template <class Tsoft>
83  void shiftToCenterOfMassFrame(Tsoft* _tsys, const PS::S64 _n) {
84  if (!pcm.is_center_shift_flag) {
85  for (int i=0; i<_n; i++) {
86  _tsys[i].pos -= pcm.pos;
87  _tsys[i].vel -= pcm.vel;
88  }
89  }
90  pcm.is_center_shift_flag = true;
91  }
92 
94 
98  template <class Tsoft>
99  void shiftToOriginFrame(Tsoft* _tsys, const PS::S64 _n) {
100  if (pcm.is_center_shift_flag) {
101  for (int i=0; i<_n; i++) {
102  _tsys[i].pos += pcm.pos;
103  _tsys[i].vel += pcm.vel;
104  }
105  }
106  pcm.is_center_shift_flag = false;
107  }
108 
110 
115  template <class Tsoft>
116  void calcCenterOfMass(Tsoft* _tsys, const PS::S64 _n, int _mode=3) {
117  PS::F64 mass = 0.0;
118  PS::F64vec pos_cm = PS::F64vec(0.0);
119  PS::F64vec vel_cm = PS::F64vec(0.0);
120 
121  if (_mode==1) { // center of the mass
122 //#pragma omp declare reduction(+:PS::F64vec:omp_out += omp_in) initializer (omp_priv=PS::F64vec(0.0))
123 //#pragma omp parallel for reduction(+:mass,pos_cm,vel_cm)
124  for (int i=0; i<_n; i++) {
125  auto& pi = _tsys[i];
126  PS::F64 mi = pi.mass;
127  mass += mi;
128 #ifdef NAN_CHECK_DEBUG
129  assert(!std::isnan(pi.vel.x));
130  assert(!std::isnan(pi.vel.y));
131  assert(!std::isnan(pi.vel.z));
132 #endif
133  pos_cm += mi*pi.pos;
134  vel_cm += mi*pi.vel;
135  }
136 
137 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
138  pcm.mass = PS::Comm::getSum(mass);
139  pcm.pos = PS::Comm::getSum(pos_cm);
140  pcm.vel = PS::Comm::getSum(vel_cm);
141 #else
142  pcm.mass = mass;
143  pcm.pos = pos_cm;
144  pcm.vel = vel_cm;
145 #endif
146  if (pcm.mass>0) {
147  pcm.pos /= pcm.mass;
148  pcm.vel /= pcm.mass;
149  }
150  }
151  else if (_mode==2) { // no mass weighted center
152 //#pragma omp declare reduction(+:PS::F64vec:omp_out += omp_in) initializer (omp_priv=PS::F64vec(0.0))
153 //#pragma omp parallel for reduction(+:mass,pos_cm,vel_cm)
154  for (int i=0; i<_n; i++) {
155  auto& pi = _tsys[i];
156  PS::F64 mi = pi.mass;
157  mass += mi;
158 #ifdef NAN_CHECK_DEBUG
159  assert(!std::isnan(pi.vel.x));
160  assert(!std::isnan(pi.vel.y));
161  assert(!std::isnan(pi.vel.z));
162 #endif
163  pos_cm += pi.pos;
164  vel_cm += pi.vel;
165  }
166 
167 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
168  pcm.mass = PS::Comm::getSum(mass);
169  pcm.pos = PS::Comm::getSum(pos_cm);
170  pcm.vel = PS::Comm::getSum(vel_cm);
171  PS::S64 n_glb = PS::Comm::getSum(_n);
172 #else
173  pcm.mass = mass;
174  pcm.pos = pos_cm;
175  pcm.vel = vel_cm;
176  PS::S64 n_glb = _n;
177 #endif
178  if (n_glb>0) {
179  pcm.pos /= PS::F64(n_glb);
180  pcm.vel /= PS::F64(n_glb);
181  }
182  }
183  else if (_mode==3) { // soft potential
184  PS::F64 pot_tot = 0.0;
185  for (int i=0; i<_n; i++) {
186  auto& pi = _tsys[i];
187  PS::F64 mi = pi.mass;
188  PS::F64 poti = pi.pot_soft;
189 #ifdef EXTERNAL_POT_IN_PTCL
190  poti -= pi.pot_ext; // remove external potential
191 #endif
192  mass += mi;
193 #ifdef NAN_CHECK_DEBUG
194  assert(!std::isnan(pi.vel.x));
195  assert(!std::isnan(pi.vel.y));
196  assert(!std::isnan(pi.vel.z));
197 #endif
198  pos_cm += poti*pi.pos;
199  vel_cm += poti*pi.vel;
200  pot_tot += poti;
201  }
202 
203 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
204  pcm.mass = PS::Comm::getSum(mass);
205  pcm.pos = PS::Comm::getSum(pos_cm);
206  pcm.vel = PS::Comm::getSum(vel_cm);
207  PS::F64 pot_tot_glb = PS::Comm::getSum(pot_tot);
208 #else
209  pcm.mass = mass;
210  pcm.pos = pos_cm;
211  pcm.vel = vel_cm;
212  PS::F64 pot_tot_glb = pot_tot;
213 #endif
214  if (pot_tot_glb!=0) {
215  pcm.pos /= pot_tot_glb;
216  pcm.vel /= pot_tot_glb;
217  }
218  }
219  }
220 
222 
227  template <class Tsoft>
228  void calcAndShiftCenterOfMass(Tsoft* _tsys, const PS::S64 _n, const int _mode=3, const bool initial_flag=false) {
229  if (initial_flag) {
230  if(pcm.is_center_shift_flag) {
231  std::cerr<<"Error: particle system is in c.m. frame, cannot initial c.m."<<std::endl;
232  abort();
233  }
234  pcm.clear();
235  pcm.is_center_shift_flag=true;
236  }
237  assert(pcm.is_center_shift_flag);
238 
239  ParticleCM pcm_bk = pcm;
240  calcCenterOfMass(_tsys, _n, _mode);
241 
242  // correct particle
243  for (int i=0; i<_n; i++) {
244  _tsys[i].pos -= pcm.pos;
245  _tsys[i].vel -= pcm.vel;
246  }
247  pcm.pos += pcm_bk.pos;
248  pcm.vel += pcm_bk.vel;
249  }
250 
252 
256  void printColumnTitle(std::ofstream & _fout, const PS::S32 _width=20) {
257  _fout<<std::setw(_width)<<"Time"
258  <<std::setw(_width)<<"N_real_loc"
259  <<std::setw(_width)<<"N_real_glb"
260  <<std::setw(_width)<<"N_all_loc"
261  <<std::setw(_width)<<"N_all_glb"
262  <<std::setw(_width)<<"N_rm_glb"
263  <<std::setw(_width)<<"N_esc_glb";
264  energy.printColumnTitle(_fout, _width);
265  pcm.printColumnTitle(_fout, _width);
266  }
267 
269 
273  void printColumn(std::ofstream & _fout, const PS::S32 _width=20) {
274  _fout<<std::setw(_width)<<time
275  <<std::setw(_width)<<n_real_loc
276  <<std::setw(_width)<<n_real_glb
277  <<std::setw(_width)<<n_all_loc
278  <<std::setw(_width)<<n_all_glb
279  <<std::setw(_width)<<n_remove_glb
280  <<std::setw(_width)<<n_escape_glb;
281  energy.printColumn(_fout, _width);
282  pcm.printColumn(_fout, _width);
283  }
284 
286 
290  void print(std::ostream & _fout, const PS::S32 _precision=7) {
291  _fout<<"Time: "<<std::setprecision(15)<<time
292  <<std::setprecision(_precision);
293  _fout<<" N_real(loc): "<<n_real_loc
294  <<" N_real(glb): "<<n_real_glb
295  <<" N_all(loc): "<<n_all_loc
296  <<" N_all(glb): "<<n_all_glb
297  <<" N_remove(glb): "<<n_remove_glb
298  <<" N_escape(glb): "<<n_escape_glb
299  <<std::endl;
300  energy.print(_fout);
301  pcm.print(_fout);
302  _fout<<std::endl;
303  }
304 };
Status::ParticleCM::vel
PS::F64vec vel
Definition: status.hpp:21
EnergyAndMomemtum::printColumn
void printColumn(std::ofstream &_fout, const PS::S32 _width=16) const
print data of class members using column style
Definition: energy.hpp:145
Status::pcm
struct Status::ParticleCM pcm
Status::ParticleCM::clear
void clear()
Definition: status.hpp:66
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
Status::n_escape_glb
PS::S32 n_escape_glb
Definition: status.hpp:15
PIKG::F64vec
Vector3< F64 > F64vec
Definition: pikg_vector.hpp:167
EnergyAndMomemtum
class for collecting and calculating the energy and angular momemtum of the system
Definition: energy.hpp:4
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
Status::ParticleCM::printColumnTitle
void printColumnTitle(std::ofstream &_fout, const PS::S32 _width=20)
print titles of class members using column style
Definition: status.hpp:31
Status::time
PS::F64 time
Definition: status.hpp:9
Status
class for measure the status of the system
Definition: status.hpp:7
Status::Status
Status()
Definition: status.hpp:75
PIKG::S64
int64_t S64
Definition: pikg_vector.hpp:23
Status::n_remove_glb
PS::S32 n_remove_glb
Definition: status.hpp:14
Status::printColumnTitle
void printColumnTitle(std::ofstream &_fout, const PS::S32 _width=20)
print titles of class members using column style
Definition: status.hpp:256
Status::n_real_glb
PS::S64 n_real_glb
Definition: status.hpp:11
Status::n_all_glb
PS::S64 n_all_glb
Definition: status.hpp:13
Status::ParticleCM::printColumn
void printColumn(std::ofstream &_fout, const PS::S32 _width=20)
print data of class members using column style
Definition: status.hpp:46
Status::ParticleCM::print
void print(std::ostream &_fout)
print title and values in one lines
Definition: status.hpp:60
Status::n_all_loc
PS::S64 n_all_loc
Definition: status.hpp:12
Status::half_mass_radius
PS::F64 half_mass_radius
Definition: status.hpp:16
Status::print
void print(std::ostream &_fout, const PS::S32 _precision=7)
print title and values in one lines
Definition: status.hpp:290
Status::ParticleCM::pos
PS::F64vec pos
Definition: status.hpp:20
Status::calcAndShiftCenterOfMass
void calcAndShiftCenterOfMass(Tsoft *_tsys, const PS::S64 _n, const int _mode=3, const bool initial_flag=false)
calculate the center of system and shift particle systems to center frame
Definition: status.hpp:228
EnergyAndMomemtum::print
void print(std::ostream &_fout=std::cout, const PS::S32 _width=16)
print title and values in one lines
Definition: energy.hpp:52
EnergyAndMomemtum::printColumnTitle
void printColumnTitle(std::ofstream &_fout, const PS::S32 _width=16) const
print titles of class members using column style
Definition: energy.hpp:111
Status::n_real_loc
PS::S64 n_real_loc
Definition: status.hpp:10
Status::energy
EnergyAndMomemtum energy
Definition: status.hpp:17
Status::shiftToOriginFrame
void shiftToOriginFrame(Tsoft *_tsys, const PS::S64 _n)
shift particle system center back to original frame
Definition: status.hpp:99
Status::ParticleCM
Definition: status.hpp:18
Status::ParticleCM::is_center_shift_flag
bool is_center_shift_flag
Definition: status.hpp:22
energy.hpp
Status::calcCenterOfMass
void calcCenterOfMass(Tsoft *_tsys, const PS::S64 _n, int _mode=3)
calculate the center of system
Definition: status.hpp:116
Status::ParticleCM::ParticleCM
ParticleCM()
Definition: status.hpp:24
Status::ParticleCM::mass
PS::F64 mass
Definition: status.hpp:19
Status::shiftToCenterOfMassFrame
void shiftToCenterOfMassFrame(Tsoft *_tsys, const PS::S64 _n)
shift particle system center to c.m. frame
Definition: status.hpp:83
Status::printColumn
void printColumn(std::ofstream &_fout, const PS::S32 _width=20)
print data of class members using column style
Definition: status.hpp:273
particle_base.hpp