PeTar
N-body code for collisional gravitational systems
energy.hpp
Go to the documentation of this file.
1 #pragma once
2 
4 class EnergyAndMomemtum{
5 public:
6  PS::F64 error_cum_pre; // previous cumulative error
7  PS::F64 ekin;
8  PS::F64 epot;
9  PS::F64 etot_ref; // energy reference (initial value) for calculating energy error
10 #ifdef HARD_CHECK_ENERGY
11  PS::F64 de_change_cum; // cumulative energy change
12  PS::F64 de_change_binary_interrupt; // cumulative energy change due to interruption (only hard, without cluster cm energy change)
13  PS::F64 de_change_modify_single; // cumulative energy change due to modification of particle (only hard, without cluster cm energy change)
14  PS::F64 error_hard_cum_pre;
15  PS::F64 error_hard_cum;
16  // slowdown energy
17  PS::F64 error_sd_cum_pre; // previous cumulative slowdown error
18  PS::F64 ekin_sd;
19  PS::F64 epot_sd;
20  PS::F64 etot_sd_ref; // slowdown energy reference (initial value) for calculating energy error
21  PS::F64 de_sd_change_cum; // cumulative slowdown energy change
22  PS::F64 de_sd_change_binary_interrupt; // cumulative slowdown energy change du to interruption (only hard, without cluster cm energy change)
23  PS::F64 de_sd_change_modify_single; // cumulative slowdown energy change due to modification of particle (only hard, without cluster cm energy change)
24  PS::F64 error_hard_sd_cum_pre;
25  PS::F64 error_hard_sd_cum;
26 #endif
27  PS::F64 error_Lt_cum_pre; // previous angular momentum error
28  PS::F64 Lt; // total angular momemtum
29  PS::F64vec L; // angular mommentum;
30  PS::F64vec L_ref; // total angular momemtum reference
31 
32  EnergyAndMomemtum() {
33  clear();
34  }
35 
36  void clear(){
37  error_cum_pre = ekin = epot = etot_ref = 0.0;
38 #ifdef HARD_CHECK_ENERGY
39  de_change_cum = de_change_binary_interrupt = de_change_modify_single = 0.0;
40  error_hard_cum = error_hard_cum_pre = 0.0;
41  error_sd_cum_pre = ekin_sd = epot_sd = etot_sd_ref = de_sd_change_cum = de_sd_change_binary_interrupt = de_sd_change_modify_single = 0.0;
42  error_hard_sd_cum = error_hard_sd_cum_pre = 0.0;
43 #endif
44  error_Lt_cum_pre = Lt = 0.0;
45  L = L_ref = PS::F64vec(0.0);
46  }
47 
49 
52  void print(std::ostream & _fout=std::cout, const PS::S32 _width=16) {
53  _fout<<"Energy: "
54  <<std::setw(_width)<<"Error/Total"
55  <<std::setw(_width)<<"Error"
56  <<std::setw(_width)<<"Error_cum"
57  <<std::setw(_width)<<"Total"
58  <<std::setw(_width)<<"Kinetic"
59  <<std::setw(_width)<<"Potential"
60 #ifdef HARD_CHECK_ENERGY
61  <<std::setw(_width)<<"Modify"
62  <<std::setw(_width)<<"Modify_group"
63  <<std::setw(_width)<<"Modify_single"
64  <<std::setw(_width)<<"Error_PP"
65  <<std::setw(_width)<<"Error_PP_cum"
66 #endif
67  <<std::endl;
68  _fout<<"Physic: "
69  <<std::setw(_width)<<(getEnergyError() - error_cum_pre)/(ekin+epot)
70  <<std::setw(_width)<<getEnergyError() - error_cum_pre
71  <<std::setw(_width)<<getEnergyError()
72  <<std::setw(_width)<<ekin + epot
73  <<std::setw(_width)<<ekin
74  <<std::setw(_width)<<epot
75 #ifdef HARD_CHECK_ENERGY
76  <<std::setw(_width)<<de_change_cum
77  <<std::setw(_width)<<de_change_binary_interrupt
78  <<std::setw(_width)<<de_change_modify_single
79  <<std::setw(_width)<<error_hard_cum - error_hard_cum_pre
80  <<std::setw(_width)<<error_hard_cum
81 #endif
82  <<std::endl;
83 #ifdef HARD_CHECK_ENERGY
84  _fout<<"Slowdown:"
85  <<std::setw(_width)<<(getEnergyErrorSlowDown() - error_sd_cum_pre)/(ekin_sd+epot_sd)
86  <<std::setw(_width)<<getEnergyErrorSlowDown() - error_sd_cum_pre
87  <<std::setw(_width)<<getEnergyErrorSlowDown()
88  <<std::setw(_width)<<ekin_sd + epot_sd
89  <<std::setw(_width)<<ekin_sd
90  <<std::setw(_width)<<epot_sd
91  <<std::setw(_width)<<de_sd_change_cum
92  <<std::setw(_width)<<de_sd_change_binary_interrupt
93  <<std::setw(_width)<<de_sd_change_modify_single
94  <<std::setw(_width)<<error_hard_sd_cum - error_hard_sd_cum_pre
95  <<std::setw(_width)<<error_hard_sd_cum
96  <<std::endl;
97 #endif
98  _fout<<"Angular Momemtum:"
99  <<" |L|err: "<<getMomentumError() - error_Lt_cum_pre
100  <<" |L|err_cum: "<<getMomentumError()
101  <<" L: "<<L
102  <<" |L|: "<<Lt
103  <<std::endl;
104  }
105 
107 
111  void printColumnTitle(std::ofstream & _fout, const PS::S32 _width=16) const {
112  _fout<<std::setw(_width)<<"Error"
113  <<std::setw(_width)<<"Error_cum"
114  <<std::setw(_width)<<"Ekin"
115  <<std::setw(_width)<<"Epot"
116  <<std::setw(_width)<<"Etot"
117 #ifdef HARD_CHECK_ENERGY
118  <<std::setw(_width)<<"dE_modify"
119  <<std::setw(_width)<<"dE_interrupt"
120  <<std::setw(_width)<<"Error_hard"
121  <<std::setw(_width)<<"Error_hard_cum"
122  <<std::setw(_width)<<"Error_SD"
123  <<std::setw(_width)<<"Error_SD_cum"
124  <<std::setw(_width)<<"Ekin_SD"
125  <<std::setw(_width)<<"Epot_SD"
126  <<std::setw(_width)<<"Etot_SD"
127  <<std::setw(_width)<<"dE_modify_SD"
128  <<std::setw(_width)<<"dE_interrupt_SD"
129  <<std::setw(_width)<<"Error_hard_SD"
130  <<std::setw(_width)<<"Error_hard_SD_cum"
131 #endif
132  <<std::setw(_width)<<"|L|error"
133  <<std::setw(_width)<<"|L|error_cum"
134  <<std::setw(_width)<<"Lx"
135  <<std::setw(_width)<<"Ly"
136  <<std::setw(_width)<<"Lz"
137  <<std::setw(_width)<<"|L|";
138  }
139 
141 
145  void printColumn(std::ofstream & _fout, const PS::S32 _width=16) const {
146  _fout<<std::setw(_width)<<getEnergyError() - error_cum_pre
147  <<std::setw(_width)<<getEnergyError()
148  <<std::setw(_width)<<ekin
149  <<std::setw(_width)<<epot
150  <<std::setw(_width)<<ekin+epot
151 #ifdef HARD_CHECK_ENERGY
152  <<std::setw(_width)<<de_change_cum
153  <<std::setw(_width)<<de_change_binary_interrupt
154  <<std::setw(_width)<<error_hard_cum - error_hard_cum_pre
155  <<std::setw(_width)<<error_hard_cum
156  <<std::setw(_width)<<getEnergyErrorSlowDown() - error_sd_cum_pre
157  <<std::setw(_width)<<getEnergyErrorSlowDown()
158  <<std::setw(_width)<<ekin_sd
159  <<std::setw(_width)<<epot_sd
160  <<std::setw(_width)<<ekin_sd+epot_sd
161  <<std::setw(_width)<<de_sd_change_cum
162  <<std::setw(_width)<<de_sd_change_binary_interrupt
163  <<std::setw(_width)<<error_hard_sd_cum - error_hard_sd_cum_pre
164  <<std::setw(_width)<<error_hard_sd_cum
165 #endif
166  <<std::setw(_width)<<getMomentumError() - error_Lt_cum_pre
167  <<std::setw(_width)<<getMomentumError()
168  <<std::setw(_width)<<L.x
169  <<std::setw(_width)<<L.y
170  <<std::setw(_width)<<L.z
171  <<std::setw(_width)<<Lt;
172  }
173 
174 #ifdef HARD_CHECK_ENERGY
175  void writeAscii(FILE* _fout) {
176  fprintf(_fout, "%26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e ",
177  error_cum_pre, ekin, epot, etot_ref, de_change_cum, de_change_binary_interrupt,
178  error_sd_cum_pre, ekin_sd, epot_sd, etot_sd_ref, de_sd_change_cum, de_sd_change_binary_interrupt,
179  error_Lt_cum_pre, Lt, L[0], L[1], L[2],
180  L_ref[0], L_ref[1], L_ref[2]);
181  }
182 #else
183  void writeAscii(FILE* _fout) {
184  fprintf(_fout, "%26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e ",
185  error_cum_pre, ekin, epot, etot_ref,
186  error_Lt_cum_pre, Lt, L[0], L[1], L[2],
187  L_ref[0], L_ref[1], L_ref[2]);
188  }
189 #endif
190 
191  void writeBinary(FILE* _fout) {
192  fwrite(&ekin, sizeof(EnergyAndMomemtum), 1, _fout);
193  }
194 
196 
202  template<class Tptcl>
203  void calc(const Tptcl* _particles,
204  const PS::S32 _n_particle,
205  const bool _init_flag=false, const PS::F64vec* _pos_offset=NULL, const PS::F64vec* _vel_offset=NULL) {
206  assert(Ptcl::group_data_mode == GroupDataMode::artificial);
207  ekin = epot = 0.0;
208  L = PS::F64vec(0.0);
209 //#pragma omp declare reduction(+:PS::F64vec:omp_out += omp_in) initializer (omp_priv=PS::F64vec(0.0))
210 //#pragma omp parallel for reduction(+:epot,ekin,L)
211  for(PS::S32 i=0; i<_n_particle; i++){
212  PS::F64 mi = _particles[i].mass;
213  auto pi_artificial = _particles[i].group_data.artificial;
214  if(pi_artificial.isMember()) mi = pi_artificial.getMassBackup();
215 #ifdef HARD_DEBUG
216  assert(_particles[i].id>0&&(pi_artificial.isMember()||pi_artificial.isSingle()));
217  assert(mi>0);
218 #endif
219 
220  PS::F64vec pi = _particles[i].pos;
221  if (_pos_offset!=NULL) pi += *_pos_offset;
222 
223  PS::F64vec vi = _particles[i].vel;
224  if (_vel_offset!=NULL) vi += *_vel_offset;
225 
226 #ifdef EXTERNAL_POT_IN_PTCL
227  epot += 0.5 * mi * (_particles[i].pot_tot + _particles[i].pot_ext);
228 #else
229  epot += 0.5 * mi * _particles[i].pot_tot;
230 #endif
231  ekin += 0.5 * mi * vi * vi;
232  L += pi ^ (mi*vi);
233  }
234  Lt = std::sqrt(L*L);
235  if (_init_flag) {
236  etot_ref = ekin + epot;
237 #ifdef HARD_CHECK_ENERGY
238  etot_sd_ref = etot_ref;
239 #endif
240  L_ref = L;
241  }
242  }
243 
245 
251  template<class Tptcl>
252  void calc(const Tptcl* _particles,
253  const PS::S32* _particle_index,
254  const PS::S32 _n_particle,
255  const bool _init_flag=false) {
256  ekin = epot = 0.0;
257  L = PS::F64vec(0.0);
258 //#pragma omp declare reduction(+:PS::F64vec:omp_out += omp_in) initializer (omp_priv=PS::F64vec(0.0))
259 //#pragma omp parallel for reduction(+:epot,ekin,L)
260  for(PS::S32 k=0; k<_n_particle; k++){
261  PS::S32 i = _particle_index[k];
262  PS::F64 mi = _particles[i].mass;
263  PS::F64vec vi = _particles[i].vel;
264 #ifdef EXTERNAL_POT_IN_PTCL
265  epot += 0.5 * mi * (_particles[i].pot_tot + _particles[i].pot_ext);
266 #else
267  epot += 0.5 * mi * _particles[i].pot_tot;
268 #endif
269  ekin += 0.5 * mi * vi * vi;
270  L += _particles[i].pos ^ (mi*vi);
271  }
272  Lt = std::sqrt(L*L);
273  if (_init_flag) {
274  etot_ref = ekin + epot;
275 #ifdef HARD_CHECK_ENERGY
276  etot_sd_ref = etot_ref;
277 #endif
278  L_ref = L;
279  }
280  }
281 
283 
286  void getSumMultiNodes(const bool _init_flag=false) {
287  ekin = PS::Comm::getSum(ekin);
288  epot = PS::Comm::getSum(epot);
289  L = PS::Comm::getSum(L);
290  Lt = std::sqrt(L*L);
291  if (_init_flag) {
292  etot_ref = ekin + epot;
293 #ifdef HARD_CHECK_ENERGY
294  etot_sd_ref = etot_ref;
295 #endif
296  L_ref = L;
297  }
298  }
299 
301  void saveEnergyError() {
302  error_cum_pre = getEnergyError();
303 #ifdef HARD_CHECK_ENERGY
304  error_sd_cum_pre = getEnergyErrorSlowDown();
305  error_hard_cum_pre = error_hard_cum;
306  error_hard_sd_cum_pre = error_hard_sd_cum;
307 #endif
308  error_Lt_cum_pre = getMomentumError();
309  }
310 
312  PS::F64 getEnergyError() const {
313  return ekin + epot - etot_ref;
314  }
315 
316 #ifdef HARD_CHECK_ENERGY
317  PS::F64 getEnergyErrorSlowDown() const {
319  return ekin_sd + epot_sd - etot_sd_ref;
320  }
321 #endif
322 
324  PS::F64 getMomentumError() const {
325  PS::F64vec dL= L- L_ref;
326  return std::sqrt(dL*dL);
327  }
328 
329  /*
330  EnergyAndMomemtum operator -(const EnergyAndMomemtum& eng){
331  EnergyAndMomemtum diff;
332  diff.ekin = ekin - eng.ekin;
333  diff.epot = epot - eng.epot;
334  diff.etot = etot - eng.etot;
335  diff.L = L - eng.L;
336  diff.Lt = std::sqrt(diff.L*diff.L);
337  return diff;
338  }
339 
340  void relative(const EnergyAndMomemtum& ref) {
341  ekin /= ref.ekin;
342  epot /= ref.epot;
343  etot /= ref.etot;
344  L /= ref.Lt;
345  Lt /= ref.Lt;
346  }
347  */
348 };
349 
EnergyAndMomemtum::ekin
PS::F64 ekin
Definition: energy.hpp:7
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
EnergyAndMomemtum::Lt
PS::F64 Lt
Definition: energy.hpp:28
EnergyAndMomemtum::clear
void clear()
Definition: energy.hpp:36
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
EnergyAndMomemtum::saveEnergyError
void saveEnergyError()
save current energy error
Definition: energy.hpp:301
EnergyAndMomemtum::L_ref
PS::F64vec L_ref
Definition: energy.hpp:30
Ptcl::group_data_mode
static GroupDataMode group_data_mode
Definition: ptcl.hpp:47
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
EnergyAndMomemtum::calc
void calc(const Tptcl *_particles, const PS::S32 *_particle_index, const PS::S32 _n_particle, const bool _init_flag=false)
calculate the system kinetic and potential energy of particles
Definition: energy.hpp:252
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
GroupDataMode::artificial
@ artificial
EnergyAndMomemtum::writeAscii
void writeAscii(FILE *_fout)
Definition: energy.hpp:183
EnergyAndMomemtum::calc
void calc(const Tptcl *_particles, const PS::S32 _n_particle, const bool _init_flag=false, const PS::F64vec *_pos_offset=NULL, const PS::F64vec *_vel_offset=NULL)
calculate the system kinetic and potential energy of particles
Definition: energy.hpp:203
EnergyAndMomemtum::error_cum_pre
PS::F64 error_cum_pre
Definition: energy.hpp:6
EnergyAndMomemtum::etot_ref
PS::F64 etot_ref
Definition: energy.hpp:9
EnergyAndMomemtum::getMomentumError
PS::F64 getMomentumError() const
get angular momemtum (value) error
Definition: energy.hpp:324
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::error_Lt_cum_pre
PS::F64 error_Lt_cum_pre
Definition: energy.hpp:27
EnergyAndMomemtum::L
PS::F64vec L
Definition: energy.hpp:29
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
EnergyAndMomemtum::getSumMultiNodes
void getSumMultiNodes(const bool _init_flag=false)
get summation of kinetic, potential energy and angular momemtum of all MPI processes
Definition: energy.hpp:286
EnergyAndMomemtum::writeBinary
void writeBinary(FILE *_fout)
Definition: energy.hpp:191
EnergyAndMomemtum::epot
PS::F64 epot
Definition: energy.hpp:8
EnergyAndMomemtum::getEnergyError
PS::F64 getEnergyError() const
get energy error
Definition: energy.hpp:312
EnergyAndMomemtum::EnergyAndMomemtum
EnergyAndMomemtum()
Definition: energy.hpp:32