SlowDown Algorithmic Regularization (SDAR)
Algorithmic Regularization with slowdown method for integrating few-body motions
information.h
Go to the documentation of this file.
1 #pragma once
2 #include "Common/Float.h"
3 #include "Common/binary_tree.h"
4 #include "AR/slow_down.h"
5 
6 namespace AR {
7 
9  class BinarySlowDown: public COMM::Binary {
10  public:
11  SlowDown slowdown;
12  Float stab_check_time;
13 
15 
17  void writeBinary(FILE *_fp) const {
18  fwrite(this, sizeof(*this),1,_fp);
19  }
20 
22 
24  void readBinary(FILE *_fin) {
25  size_t rcount = fread(this, sizeof(*this),1,_fin);
26  if (rcount<1) {
27  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
28  abort();
29  }
30  }
31 
33 
37  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
38  Binary::printColumnTitle(_fout, _width);
39  SlowDown::printColumnTitle(_fout,_width);
40  }
41 
43 
47  void printColumn(std::ostream & _fout, const int _width=20){
48  Binary::printColumn(_fout, _width);
49  slowdown.printColumn(_fout,_width);
50  }
51 
53 
55  //void writeAscii(std::ostream& _fout) const {
56  // COMM::BinaryTree<Tparticle>::writeAscii(_fout);
57  // SlowDown::writeAscii(_fout);
58  //}
59 
61 
63  //void readAscii(std::istream& _fin) {
64  // COMM::BinaryTree<Tparticle>::readAscii(_fin);
65  // SlowDown::readAscii(_fin);
66  //}
67 
68  };
69 
71  template <class Tparticle>
72  using BinaryTree=COMM::BinaryTree<Tparticle,BinarySlowDown>;
73 
75 
79  enum class FixStepOption {always, later, none};
80 
82 
89  template <class Tparticle, class Tpcm>
90  class Information{
91  public:
92  Float ds;
93  Float time_offset;
94  Float r_break_crit; // group break radius criterion
95  FixStepOption fix_step_option;
96  COMM::List<BinaryTree<Tparticle>> binarytree;
97 #ifdef AR_DEBUG_DUMP
98  bool dump_flag;
99 #endif
100 
102  Information(): ds(0.0), time_offset(0.0), r_break_crit(-1.0), fix_step_option(AR::FixStepOption::none), binarytree() {
103 #ifdef AR_DEBUG_DUMP
104  dump_flag = false;
105 #endif
106  }
107 
109 
111  bool checkParams() {
112  ASSERT(r_break_crit>=0.0);
113  ASSERT(binarytree.getSize()>0);
114  return true;
115  }
116 
118  void reserveMem(const int _nmax) {
119  binarytree.setMode(COMM::ListMode::local);
120  binarytree.reserveMem(_nmax);
121  }
122 
124  BinaryTree<Tparticle>& getBinaryTreeRoot() const {
125  int n = binarytree.getSize();
126  ASSERT(n>0);
127  return binarytree[n-1];
128  }
129 
130  inline Float calcDsElliptic(BinaryTree<Tparticle>& _bin, const Float& _G) {
131  //kepler orbit, step ds=dt*G*m1*m2/r estimation (1/32 orbit): 2*pi/32*sqrt(G*semi/(m1+m2))*m1*m2
132  return 0.19634954084*sqrt(_G*_bin.semi/(_bin.m1+_bin.m2))*(_bin.m1*_bin.m2);
133  }
134 
135  inline Float calcDsHyperbolic(BinaryTree<Tparticle>& _bin, const Float& _G) {
136  //hyperbolic orbit, step ds=dt*m1*m2/r estimation (1/256 orbit): pi/128*sqrt(G*semi/(m1+m2))*m1*m2
137  return 0.0245436926*sqrt(-_G*_bin.semi/(_bin.m1+_bin.m2))*(_bin.m1*_bin.m2);
138  }
139 
141  void calcDsMinKeplerIter(Float& _ds_over_ebin_min_bin, Float& _ds_min_bin, Float& _ds_min_hyp, Float& _etot_sd, const Float& _G, const Float& _nest_sd_up, BinaryTree<Tparticle>& _bin, const int _intergrator_order) {
142  Float nest_sd = _nest_sd_up * _bin.slowdown.getSlowDownFactor();
143  // perturbation ratio
144  Float pert_ratio = (_bin.slowdown.pert_out>0&&_bin.slowdown.pert_in>0)? _bin.slowdown.pert_in/_bin.slowdown.pert_out: 1.0;
145  // scale step based on perturbation and sym method order
146  Float scale_factor = std::min(Float(1.0),pow(1e-1*pert_ratio,1.0/Float(_intergrator_order)));
147  for (int k=0; k<2; k++) {
148  if (_bin.isMemberTree(k)) {
149  calcDsMinKeplerIter(_ds_over_ebin_min_bin, _ds_min_bin, _ds_min_hyp, _etot_sd, _G, nest_sd, *_bin.getMemberAsTree(k), _intergrator_order);
150  }
151  }
152  // zero mass cause ds=0
153  if (_bin.m1>0&&_bin.m2>0) {
154  if (_bin.semi>0) {
155  Float dsi = calcDsElliptic(_bin, _G);
156  // scale by /Ebin_sd
157  Float ebin_sd = _G*(_bin.m1*_bin.m2)/(2*_bin.semi*nest_sd);
158  ASSERT(dsi>0&&ebin_sd>0);
159  Float ds_over_ebin = dsi*scale_factor/ebin_sd;
160  if (ds_over_ebin<_ds_over_ebin_min_bin) {
161  _ds_over_ebin_min_bin = ds_over_ebin;
162  _ds_min_bin = dsi*scale_factor;
163  }
164  _etot_sd += ebin_sd;
165  }
166  else {
167  Float dsi = calcDsHyperbolic(_bin, _G);
168  ASSERT(dsi>0);
169  //Float factor = std::min(Float(1.0), pow(nest_sd_org,Float(1.0/3.0)));
170  _ds_min_hyp = std::min(dsi, _ds_min_hyp);
171  }
172  }
173  }
174 
176 
178  Float calcDsKeplerBinaryTree(BinaryTree<Tparticle>& _bin, const int _int_order, const Float& _G, const Float& _ds_scale) {
179  Float ds_over_ebin_min=NUMERIC_FLOAT_MAX;
180  Float ds_min_hyp=NUMERIC_FLOAT_MAX;
181  Float ds_min_bin=NUMERIC_FLOAT_MAX;
182  Float etot_sd = 0.0;
183  calcDsMinKeplerIter(ds_over_ebin_min, ds_min_bin, ds_min_hyp, etot_sd, _G, 1.0, _bin, _int_order);
184  //Float ds_min_bin = ds_over_ebin_min*etot_sd/(bin_root.getMemberN()-1);
185  ASSERT(ds_min_hyp<NUMERIC_FLOAT_MAX||ds_min_bin<NUMERIC_FLOAT_MAX);
186  return _ds_scale*std::min(ds_min_bin,ds_min_hyp);
187  }
188 
190 
195  void calcDsAndStepOption(const int _int_order, const Float& _G, const Float& _ds_scale) {
196  auto& bin_root = getBinaryTreeRoot();
197  ds = calcDsKeplerBinaryTree(bin_root, _int_order, _G, _ds_scale);
198  ASSERT(ds>0);
199 
200  // Avoid too small step
201  //if (_sd_org<1.0) ds *= std::max(1.0/8.0*pow(_sd_org, 1.0/Float(_int_order)),0.125);
202  //auto& bin_root = getBinaryTreeRoot();
203  const int n_particle = bin_root.getMemberN();
204 
205  // determine the fix step option
206  //fix_step_option = FixStepOption::later;
207  fix_step_option = FixStepOption::none;
209  if (n_particle==2||bin_root.stab<1.0) fix_step_option = AR::FixStepOption::later;
211  //if (n_particle>2) {
212  // Float apo_in_max = 0;
213  // for (int j=0; j<2; j++) {
214  // if (bin_root.isMemberTree(j)) {
215  // auto* bin_sub = (COMM::BinaryTree<Tparticle>*) bin_root.getMember(j);
216  // apo_in_max = std::max(apo_in_max,bin_sub->semi*(1+bin_sub->ecc));
217  // }
218  // }
219  // Float peri_out = bin_root.semi*(1-bin_root.ecc);
220  // if (peri_out>3*apo_in_max) fix_step_option = FixStepOption::later;
221  //}
222  }
223 
225 
229  void generateBinaryTree(COMM::ParticleGroup<Tparticle, Tpcm>& _particles, const Float _G) {
230  const int n_particle = _particles.getSize();
231  ASSERT(n_particle>1);
232  binarytree.resizeNoInitialize(n_particle-1);
233  int particle_index_local[n_particle];
234  int particle_index_unused[n_particle];
235  int n_particle_real = 0;
236  int n_particle_unused=0;
237  for (int i=0; i<n_particle; i++) {
238  if (_particles[i].mass>0.0) particle_index_local[n_particle_real++] = i;
239  else particle_index_unused[n_particle_unused++] = i;
240  }
241  if (n_particle_real>1)
242  BinaryTree<Tparticle>::generateBinaryTree(binarytree.getDataAddress(), particle_index_local, n_particle_real, _particles.getDataAddress(), _G);
243 
244  // Add unused particles to the outmost orbit
245  if (n_particle_real==0) {
246  binarytree[0].setMembers(&(_particles[0]), &(_particles[1]), 0, 1);
247  binarytree[0].mass = 0.0;
248  binarytree[0].m1 = 0.0;
249  binarytree[0].m2 = 0.0;
250  for (int i=2; i<n_particle; i++) {
251  binarytree[i-1].setMembers((Tparticle*)&(binarytree[i-2]), &( _particles[i]), -1, i);
252  binarytree[i-1].mass = 0.0;
253  binarytree[i-1].m1 = 0.0;
254  binarytree[i-1].m2 = 0.0;
255  }
256  }
257  else if (n_particle_real==1) {
258  int i1 = particle_index_local[0];
259  int i2 = particle_index_unused[0];
260  binarytree[0].setMembers(&(_particles[i1]), &(_particles[i2]), i1 ,i2);
261  binarytree[0].m1 = _particles[i1].mass;
262  binarytree[0].m2 = _particles[i2].mass;
263  binarytree[0].mass = binarytree[0].m1 + binarytree[0].m2;
264  binarytree[0].pos[0] = _particles[i1].pos[0];
265  binarytree[0].pos[1] = _particles[i1].pos[1];
266  binarytree[0].pos[2] = _particles[i1].pos[2];
267  binarytree[0].vel[0] = _particles[i1].vel[0];
268  binarytree[0].vel[1] = _particles[i1].vel[1];
269  binarytree[0].vel[2] = _particles[i1].vel[2];
270  for (int i=1; i<n_particle_unused; i++) {
271  int k = particle_index_unused[i];
272  binarytree[i].setMembers((Tparticle*)&(binarytree[i-1]), &(_particles[k]), -1, k);
273  binarytree[i].m1 = binarytree[i-1].mass;
274  binarytree[i].m2 = _particles[k].mass;
275  binarytree[i].mass = binarytree[i].m1 + binarytree[i].m2;
276  binarytree[i].pos[0] = binarytree[i-1].pos[0];
277  binarytree[i].pos[1] = binarytree[i-1].pos[1];
278  binarytree[i].pos[2] = binarytree[i-1].pos[2];
279  binarytree[i].vel[0] = binarytree[i-1].vel[0];
280  binarytree[i].vel[1] = binarytree[i-1].vel[1];
281  binarytree[i].vel[2] = binarytree[i-1].vel[2];
282  }
283  }
284  else {
285  for (int i=0; i<n_particle_unused; i++) {
286  int ilast = n_particle_real-1+i;
287  ASSERT(ilast<n_particle-1);
288  int k = particle_index_unused[i];
289  binarytree[ilast].setMembers((Tparticle*)&(binarytree[ilast-1]), &(_particles[k]), -1, k);
290  binarytree[ilast].m1 = binarytree[ilast-1].mass;
291  binarytree[ilast].m2 = _particles[k].mass;
292  binarytree[ilast].mass = binarytree[ilast].m1 + binarytree[ilast].m2;
293  binarytree[ilast].pos[0] = binarytree[ilast-1].pos[0];
294  binarytree[ilast].pos[1] = binarytree[ilast-1].pos[1];
295  binarytree[ilast].pos[2] = binarytree[ilast-1].pos[2];
296  binarytree[ilast].vel[0] = binarytree[ilast-1].vel[0];
297  binarytree[ilast].vel[1] = binarytree[ilast-1].vel[1];
298  binarytree[ilast].vel[2] = binarytree[ilast-1].vel[2];
299  }
300  }
301  }
302 
304 
308  bool checkAndSetBinaryPairIDIter(BinaryTree<Tparticle>& _bin, const bool _reset_flag) {
309  bool return_flag=true;
310  Tparticle* p[2] = {_bin.getLeftMember(), _bin.getRightMember()};
311 
312  for (int i=0; i<2; i++) {
313  if (_bin.isMemberTree(i)) {
314  return_flag = return_flag & checkAndSetBinaryPairIDIter(*_bin.getMemberAsTree(i),_reset_flag);
315  }
316  }
317  for (int i=0; i<2; i++) {
318  if (!_bin.isMemberTree(i)) {
319  auto pair_id = p[1-i]->id;
320  return_flag = return_flag & (p[i]->getBinaryPairID()==pair_id);
321  if (_reset_flag) p[i]->setBinaryPairID(0);
322  else p[i]->setBinaryPairID(pair_id);
323  }
324  }
325  if (p[0]->id<p[1]->id) _bin.id = p[0]->id;
326  else _bin.id = p[1]->id;
327  return return_flag;
328  }
329 
331  void clear() {
332  ds=0.0;
333  time_offset = 0.0;
334  r_break_crit=-1.0;
335  fix_step_option = FixStepOption::none;
336  binarytree.clear();
337 #ifdef AR_DEBUG_DUMP
338  dump_flag=false;
339 #endif
340  }
341 
343 
347  void printColumnTitle(std::ostream & _fout, const int _width=20) {
348  _fout<<std::setw(_width)<<"ds";
349  _fout<<std::setw(_width)<<"Time_offset";
350  _fout<<std::setw(_width)<<"r_break_crit";
351  }
352 
354 
358  void printColumn(std::ostream & _fout, const int _width=20){
359  _fout<<std::setw(_width)<<ds;
360  _fout<<std::setw(_width)<<time_offset;
361  _fout<<std::setw(_width)<<r_break_crit;
362  }
363 
365 
367  void writeBinary(FILE *_fout) const {
368  fwrite(&ds, sizeof(int),1,_fout);
369  fwrite(&time_offset, sizeof(Float),1,_fout);
370  fwrite(&r_break_crit, sizeof(Float),1,_fout);
371  fwrite(&fix_step_option, sizeof(FixStepOption),1,_fout);
372  }
373 
375 
377  void readBinary(FILE *_fin) {
378  size_t rcount = fread(&ds, sizeof(int),1,_fin);
379  if (rcount<1) {
380  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
381  abort();
382  }
383  rcount = fread(&time_offset, sizeof(Float),1,_fin);
384  if (rcount<1) {
385  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
386  abort();
387  }
388  rcount = fread(&r_break_crit, sizeof(Float),1,_fin);
389  if (rcount<1) {
390  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
391  abort();
392  }
393  rcount = fread(&fix_step_option, sizeof(FixStepOption),1,_fin);
394  if (rcount<1) {
395  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
396  abort();
397  }
398  }
399  };
400 
401 }
AR::Information::checkAndSetBinaryPairIDIter
bool checkAndSetBinaryPairIDIter(BinaryTree< Tparticle > &_bin, const bool _reset_flag)
check binary tree member pair id, if consisent, return ture. otherwise set the member pair id
Definition: information.h:308
COMM::List< Tparticle >::getDataAddress
Tparticle * getDataAddress() const
return member data array address
Definition: list.h:170
AR::FixStepOption::later
@ later
AR::Information::calcDsAndStepOption
void calcDsAndStepOption(const int _int_order, const Float &_G, const Float &_ds_scale)
calculate ds from the inner most binary with minimum period, determine the fix step option
Definition: information.h:195
AR::BinarySlowDown::readBinary
void readBinary(FILE *_fin)
read class data to file with binary format
Definition: information.h:24
AR::Information::printColumnTitle
void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: information.h:347
AR::BinarySlowDown::writeBinary
void writeBinary(FILE *_fp) const
write class data to file with binary format
Definition: information.h:17
AR::Information
A class contains information (e.g. parameters, binary tree, indices) about the particle group.
Definition: information.h:90
Float.h
AR::BinarySlowDown::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: information.h:47
AR::Information::calcDsHyperbolic
Float calcDsHyperbolic(BinaryTree< Tparticle > &_bin, const Float &_G)
Definition: information.h:135
AR
Algorithmic regularization (time transformed explicit symplectic integrator) namespace.
Definition: force.h:5
NUMERIC_FLOAT_MAX
const Float NUMERIC_FLOAT_MAX
Definition: Float.h:29
slow_down.h
AR::FixStepOption::always
@ always
COMM::Binary
Binary parameter class.
Definition: binary_tree.h:12
AR::Information::calcDsKeplerBinaryTree
Float calcDsKeplerBinaryTree(BinaryTree< Tparticle > &_bin, const int _int_order, const Float &_G, const Float &_ds_scale)
calculate average kepler ds iterately for a binary tree
Definition: information.h:178
AR::SlowDown::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: slow_down.h:193
COMM::BinaryTree::getLeftMember
Tptcl * getLeftMember() const
get left member
Definition: binary_tree.h:963
binary_tree.h
AR::Information::ds
Float ds
Definition: information.h:92
AR::BinarySlowDown::slowdown
SlowDown slowdown
Definition: information.h:11
AR::Information::Information
Information()
a list of binary tree that contain the hierarchical orbital parameters of the particle group.
Definition: information.h:102
AR::Information::binarytree
COMM::List< BinaryTree< Tparticle > > binarytree
fix step option for integration
Definition: information.h:96
AR::Information::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: information.h:358
COMM::BinaryTree::generateBinaryTree
static void generateBinaryTree(BinaryTreeLocal _bins[], int _ptcl_list[], const int _n, Tptcl *_ptcl, const Float _G)
Generate kepler binary tree for a group of particles.
Definition: binary_tree.h:666
AR::Information::generateBinaryTree
void generateBinaryTree(COMM::ParticleGroup< Tparticle, Tpcm > &_particles, const Float _G)
generate binary tree for the particle group
Definition: information.h:229
COMM::ParticleGroup
Particle group class to store and manage a group of particle.
Definition: particle_group.h:13
COMM::List< Tparticle >::getSize
int getSize() const
get current member number
Definition: list.h:151
Float
double Float
Definition: Float.h:25
AR::Information::readBinary
void readBinary(FILE *_fin)
read class data to file with binary format
Definition: information.h:377
AR::Information::calcDsMinKeplerIter
void calcDsMinKeplerIter(Float &_ds_over_ebin_min_bin, Float &_ds_min_bin, Float &_ds_min_hyp, Float &_etot_sd, const Float &_G, const Float &_nest_sd_up, BinaryTree< Tparticle > &_bin, const int _intergrator_order)
iteration function to calculate average kepler ds for a binary tree
Definition: information.h:141
COMM::BinaryTree::isMemberTree
bool isMemberTree(const size_t i) const
Definition: binary_tree.h:957
AR::BinarySlowDown::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: information.h:37
COMM::BinaryTree::getMemberAsTree
BinaryTreeLocal * getMemberAsTree(const size_t i) const
get member as tree
Definition: binary_tree.h:946
COMM::BinaryTree::getRightMember
Tptcl * getRightMember() const
get right member
Definition: binary_tree.h:973
AR::FixStepOption
FixStepOption
Fix step options for integration with adjusted step (not for time sychronizatio phase)
Definition: information.h:79
AR::Information::checkParams
bool checkParams()
check whether parameters values are correct initialized
Definition: information.h:111
AR::Information::fix_step_option
FixStepOption fix_step_option
Definition: information.h:95
AR::Information::writeBinary
void writeBinary(FILE *_fout) const
write class data to file with binary format
Definition: information.h:367
COMM::List
list class to store and manage a group of member
Definition: list.h:19
AR::Information::calcDsElliptic
Float calcDsElliptic(BinaryTree< Tparticle > &_bin, const Float &_G)
Definition: information.h:130
AR::BinarySlowDown
binary parameter with slowdown
Definition: information.h:9
AR::FixStepOption::none
@ none
COMM::BinaryTree
Binary tree cell.
Definition: binary_tree.h:492
AR::Information::clear
void clear()
clear function
Definition: information.h:331
COMM::ListMode::local
@ local
AR::Information::getBinaryTreeRoot
BinaryTree< Tparticle > & getBinaryTreeRoot() const
get the root of binary tree
Definition: information.h:124
AR::BinarySlowDown::stab_check_time
Float stab_check_time
Definition: information.h:12
AR::Information::r_break_crit
Float r_break_crit
offset of time to obtain real physical time (real time = TimeTransformedSymplecticIntegrator:time_ + ...
Definition: information.h:94
AR::SlowDown
Slow-down parameter control class.
Definition: slow_down.h:11
AR::Information::reserveMem
void reserveMem(const int _nmax)
reserve memory of binarytree list
Definition: information.h:118
AR::Information::time_offset
Float time_offset
initial step size for integration
Definition: information.h:93
AR::SlowDown::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: slow_down.h:204