PeTar
N-body code for collisional gravitational systems
petar.hpp
Go to the documentation of this file.
1 #pragma once
2 #ifdef P3T_64BIT
3 #define CALC_EP_64bit
4 #define CALC_SP_64bit
5 #define RSQRT_NR_EPJ_X4
6 #define RSQRT_NR_SPJ_X4
7 
8 #elif P3T_MIXBIT
9 #define CALC_EP_64bit
10 #define RSQRT_NR_EPJ_X4
11 
12 #else
13 #define RSQRT_NR_EPJ_X2
14 //#define RSQRT_NR_SPJ_X2
15 #endif
16 
17 #if defined(INTRINSIC_K) || defined(INTRINSIC_X86)
18 #define INTRINSIC
19 #endif
20 
21 #include<iostream>
22 #include<iomanip>
23 #include<fstream>
24 #include<string>
25 #include<sstream>
26 //#include<unistd.h>
27 #include<getopt.h>
28 
29 #ifdef MPI_DEBUG
30 #include <mpi.h>
31 // Send recv debug
32 int MPI_Isend(void* buffer, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request* req)
33 {
34  int ret;
35 
36  int size;
37  MPI_Type_size(datatype, &size);
38  std::cerr<<"MPI_ISend count "<<count<<" dest "<<dest<<" datatype "<<size<<" tag "<<tag<<std::endl;
39 
40  ret = PMPI_Isend(buffer, count, datatype, dest, tag, comm, req);
41 
42  return ret;
43 }
44 
45 // Debug
46 int MPI_Irecv(void* buffer, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request* req)
47 {
48  int ret;
49 
50  int size;
51  MPI_Type_size(datatype, &size);
52  std::cerr<<"MPI_IRecv count "<<count<<" dest "<<dest<<" datatype "<<size<<" tag "<<tag<<std::endl;
53 
54  ret = PMPI_Irecv(buffer, count, datatype, dest, tag, comm, req);
55 
56  return ret;
57 }
58 #endif
59 
60 #include<get_version.hpp>
61 #include<particle_simulator.hpp>
62 #include"hard_assert.hpp"
63 #include"soft_ptcl.hpp"
64 #include"soft_force.hpp"
65 #ifdef USE_GPU
66 #include"force_gpu_cuda.hpp"
67 #endif
68 #ifdef USE_FUGAKU
69 #include "force_fugaku.hpp"
70 #endif
71 #include"energy.hpp"
72 #include"hard.hpp"
73 #include"io.hpp"
74 #include"status.hpp"
75 #include"particle_distribution_generator.hpp"
76 #include"domain.hpp"
77 #include"cluster_list.hpp"
78 #include"kickdriftstep.hpp"
79 #ifdef PROFILE
80 #include"profile.hpp"
81 #endif
82 #include"static_variables.hpp"
83 #include"escaper.hpp"
84 #ifdef GALPY
85 #include"galpy_interface.h"
86 #endif
87 
89 class IOParamsPeTar{
90 public:
91  // IO parameters
92  IOParamsContainer input_par_store;
93  IOParams<PS::F64> ratio_r_cut;
94  IOParams<PS::F64> theta;
95  IOParams<PS::S64> n_leaf_limit;
96  IOParams<PS::S64> n_group_limit;
97  IOParams<PS::S64> n_interrupt_limit;
98  IOParams<PS::S64> n_smp_ave;
99 #ifdef ORBIT_SAMPLING
100  IOParams<PS::S64> n_split;
101 #endif
102  IOParams<PS::S64> n_bin;
103  IOParams<PS::S64> n_step_per_orbit;
104  IOParams<PS::F64> time_end;
105  IOParams<PS::F64> eta;
106  IOParams<PS::F64> gravitational_constant;
107  IOParams<PS::S64> unit_set;
108  IOParams<PS::S64> n_glb;
109  IOParams<PS::S64> id_offset;
110  IOParams<PS::F64> dt_soft;
111  IOParams<PS::F64> dt_snap;
112  IOParams<PS::F64> nstep_dt_soft_kepler;
113  IOParams<PS::F64> search_vel_factor;
114  IOParams<PS::F64> search_peri_factor;
115  IOParams<PS::F64> dt_limit_hard_factor;
116  IOParams<PS::S64> dt_min_hermite_index;
117  //IOParams<PS::S64> dt_min_ar_index;
118  //IOParams<PS::F64> dt_err_pert;
119  //IOParams<PS::F64> dt_err_soft;
120  IOParams<PS::F64> e_err_ar;
121 #ifdef HARD_CHECK_ENERGY
122  IOParams<PS::F64> e_err_hard;
123 #endif
124  IOParams<PS::S64> step_limit_ar;
125  IOParams<PS::F64> eps;
126  IOParams<PS::F64> r_out;
127  IOParams<PS::F64> r_bin;
128 // IOParams<PS::F64> r_search_max;
129  IOParams<PS::F64> r_search_min;
130  IOParams<PS::F64> r_escape;
131  IOParams<PS::F64> sd_factor;
132  IOParams<PS::S64> data_format;
133  IOParams<PS::S64> write_style;
134 #ifdef STELLAR_EVOLUTION
135  IOParams<PS::S64> stellar_evolution_option;
136 #endif
137  IOParams<PS::S64> interrupt_detection_option;
138 #ifdef ADJUST_GROUP_PRINT
139  IOParams<PS::S64> adjust_group_write_option;
140 #endif
141  IOParams<PS::S64> append_switcher;
142  IOParams<std::string> fname_snp;
143  IOParams<std::string> fname_par;
144  IOParams<std::string> fname_inp;
145 
146  // flag
147  bool print_flag;
148  bool update_changeover_flag;
149  bool update_rsearch_flag;
150 
151  IOParamsPeTar(): input_par_store(),
152  ratio_r_cut (input_par_store, 0.1, "r-ratio", "r_in / r_out"),
153  theta (input_par_store, 0.3, "T", "Particle-tree opening angle theta"),
154  n_leaf_limit (input_par_store, 20, "number-leaf-limit", "Particle-tree leaf number limit", "Optimal value should be slightly >= 11 + N_bin_sample (20)"),
155 #ifdef USE__AVX512
156  n_group_limit (input_par_store, 1024, "number-group-limit", "Particle-tree group number limit", "Optimized for x86-AVX512 (1024)"),
157 #else
158  n_group_limit (input_par_store, 512, "number-group-limit", "Particle-tree group number limit", "Optimized for x86-AVX2 (512)"),
159 #endif
160  n_interrupt_limit(input_par_store, 128, "number-interrupt-limit", "Interrupted hard integrator limit"),
161  n_smp_ave (input_par_store, 100, "number-sample-average", "Average target number of sample particles per process"),
162 #ifdef ORBIT_SAMPLING
163  n_split (input_par_store, 4, "number-split", "Number of binary sample points for tree perturbation force"),
164 #endif
165  n_bin (input_par_store, 0, "b", "Number of primordial binaries for initialization (assuming the binaries' IDs are 1,2*n_bin)"),
166  n_step_per_orbit (input_par_store, 8, "number-step-tt", "Number of steps per slow-down binary orbits (binary period/tree timestep) for isolated binaries; also the maximum criterion for activating tidal tensor method"),
167  time_end (input_par_store, 10.0, "t", "End time of simulation"),
168  eta (input_par_store, 0.1, "hermite-eta", "Hermite timestep coefficient eta"),
169  gravitational_constant(input_par_store, 1.0, "G", "Gravitational constant"),
170  unit_set (input_par_store, 0, "u", "Input data unit, 0: unknown, referring to G; 1: mass:Msun, length:pc, time:Myr, velocity:pc/Myr"),
171  n_glb (input_par_store, 100000, "n", "Total number of particles, used only for a test with the internal equal-mass Plummer model generator (assuming G=1 and the input data filename is __Plummer)"),
172  id_offset (input_par_store, -1, "id-offset", "Starting ID for artificial particles, total number of real particles must always be smaller than this","n_glb+1"),
173  dt_soft (input_par_store, 0.0, "s", "Tree timestep (dt_soft), if the value is zero (default) and --nstep-dt-soft-kepler is not used, then dt_soft = 0.1*r_out/sigma_1D"),
174  dt_snap (input_par_store, 1.0, "o", "Output time interval for particle dataset snapshots"),
175  nstep_dt_soft_kepler (input_par_store, 0.0, "nstep-dt-soft-kepler", "Determines the tree timestep by P(r_in)/nstep, where P(r_in) is the binary period with the semi-major axis of r_in, nstep is the argument of this option (e.g., 32.0)", "not used"),
176  search_vel_factor(input_par_store, 3.0, "search-vel-factor", "Neighbor search coefficient for velocity check (v*dt)"),
177  search_peri_factor (input_par_store, 1.5, "search-peri-factor", "Neighbor search coefficient for periapsis check"),
178  dt_limit_hard_factor(input_par_store, 4.0, "dt-max-factor", "Limit of tree timestep/hard timestep"),
179  dt_min_hermite_index(input_par_store, 40, "dt-min-hermite", "Power index n for the smallest timestep (0.5^n) allowed in the Hermite integrator"),
180  //dt_min_ar_index (input_par_store, 64, "dt-min-ar", "Power index n for the smallest timestep (0.5^n) allowed in the ARC integrator, suppressed"),
181  //dt_err_pert (input_par_store, 1e-6, "dt-error-pert", "Maximum time synchronization error (relative) for perturbed ARC integrator, suppressed"),
182  //dt_err_soft (input_par_store, 1e-3, "dt-error-iso", "Maximum time synchronization error (relative) for no-perturber (only soft perturbation) ARC integrator, suppressed"),
183  e_err_ar (input_par_store, 1e-8, "energy-err-ar", "Maximum energy error allowed for the ARC integrator"),
184 #ifdef HARD_CHECK_ENERGY
185  e_err_hard (input_par_store, 1e-4, "energy-err-hard", "Maximum energy error allowed for the hard integrator"),
186 #endif
187  step_limit_ar(input_par_store, 1000000, "step-limit-ar", "Maximum step allowed for the ARC sym integrator"),
188  eps (input_par_store, 0.0, "soft-eps", "Softening epsilon"),
189  r_out (input_par_store, 0.0, "r", "Outer boundary radius for the changeover function (r_out), if value is zero and -s is not used, use 0.1 GM/[N^(1/3) sigma_3D^2]; if -s is given, calculate r_out from dt_soft"),
190  r_bin (input_par_store, 0.0, "r-bin", "Tidal tensor box size and the radial criterion for detecting multiple systems (binaries, triples, etc.), if value is zero, use 0.8*r_in"),
191 // r_search_max (input_par_store, 0.0, "Maximum search radius criterion", "5*r_out"),
192  r_search_min (input_par_store, 0.0, "r-search-min", "Minimum neighbor search radius for hard clusters","auto"),
193  r_escape (input_par_store, PS::LARGE_FLOAT, "r-escape", "Escape radius criterion, 0: no escaper removal; <0: remove particles when r>-r_escape; >0: remove particles when r>r_escape and energy>0"),
194  sd_factor (input_par_store, 1e-4, "slowdown-factor", "Slowdown perturbation criterion"),
195  data_format (input_par_store, 1, "i", "Data read(r)/write(w) format BINARY(B)/ASCII(A): r-B/w-A (3), r-A/w-B (2), rw-A (1), rw-B (0)"),
196  write_style (input_par_store, 1, "w", "File writing style: 0, no output; 1. write snapshots, status, and profile separately; 2. write snapshot and status in one line per step (no MPI support); 3. write only status and profile"),
197 #ifdef STELLAR_EVOLUTION
198 #ifdef BSE_BASE
199  stellar_evolution_option (input_par_store, 1, "stellar-evolution", "Stellar evolution of stars in Hermite+SDAR: 0: off; >=1: using SSE/BSE based codes; ==2: activate dynamical tide and hyperbolic gravitational wave radiation"),
200  interrupt_detection_option(input_par_store, 1, "detect-interrupt", "Stellar evolution of binaries in SDAR: 0: off; 1: using BSE based code (if '--stellar-evolution != 0)"),
201 #else
202  stellar_evolution_option (input_par_store, 0, "stellar-evolution", "Not implemented"),
203  interrupt_detection_option(input_par_store, 0, "detect-interrupt", "Interrupt integration in SDAR: 0: turn off; 1: merge two particles if their surfaces overlap; 2. merge two particles and also interrupt the hard integration"),
204 #endif
205 #else
206  interrupt_detection_option(input_par_store, 0, "detect-interrupt", "Modify orbits of AR groups based on the interruption function: 0: turn off; 1: modify inside AR integration and accumulate energy change; 2. modify and also interrupt the hard drift"),
207 #endif
208 #ifdef ADJUST_GROUP_PRINT
209  adjust_group_write_option(input_par_store, 1, "write-group-info", "Print new and end of groups: 0: no print; 1: print to file [data filename prefix].group.[MPI rank] if -w >0"),
210 #endif
211  append_switcher(input_par_store, 1, "a", "Data output style: 0 - create new output files and overwrite existing ones except snapshots; 1 - append new data to existing files"),
212  fname_snp(input_par_store, "data", "f", "Prefix of filenames for output data: [prefix].**"),
213  fname_par(input_par_store, "input.par", "p", "Input parameter file (this option should be used first before any other options)"),
214  fname_inp(input_par_store, "__NONE__", "snap-filename", "Input data file", NULL, false),
215  print_flag(false), update_changeover_flag(false), update_rsearch_flag(false) {}
216 
217 
219 
225  int read(int argc, char *argv[], const int opt_used_pre=0) {
226  static int petar_flag=-1;
227  static struct option long_options[] = {
228 #ifdef ORBIT_SAMPLING
229  {n_split.key, required_argument, &petar_flag, 0},
230 #endif
231  {search_vel_factor.key, required_argument, &petar_flag, 1},
232  {dt_limit_hard_factor.key, required_argument, &petar_flag, 2},
233  {dt_min_hermite_index.key, required_argument, &petar_flag, 3},
234  {n_group_limit.key, required_argument, &petar_flag, 4},
235  {n_leaf_limit.key, required_argument, &petar_flag, 5},
236  {n_smp_ave.key, required_argument, &petar_flag, 6},
237  {e_err_ar.key, required_argument, &petar_flag, 7},
238  {eps.key, required_argument, &petar_flag, 8},
239  {sd_factor.key, required_argument, &petar_flag, 9},
240  {ratio_r_cut.key, required_argument, &petar_flag, 10},
241  {r_bin.key, required_argument, &petar_flag, 11},
242  {search_peri_factor.key, required_argument, &petar_flag, 12},
243  {eta.key, required_argument, &petar_flag, 13},
244 #ifdef HARD_CHECK_ENERGY
245  {e_err_hard.key, required_argument, &petar_flag, 14},
246 #endif
247  {step_limit_ar.key, required_argument, &petar_flag, 15},
248  {"disable-print-info", no_argument, &petar_flag, 16},
249  {n_interrupt_limit.key, required_argument, &petar_flag, 17},
250  {interrupt_detection_option.key, required_argument, &petar_flag, 18},
251  {n_step_per_orbit.key, required_argument, &petar_flag, 19},
252 #ifdef STELLAR_EVOLUTION
253  {stellar_evolution_option.key, required_argument, &petar_flag, 20},
254 #endif
255  {r_escape.key, required_argument, &petar_flag, 21},
256  {r_search_min.key, required_argument, &petar_flag, 22},
257  {id_offset.key, required_argument, &petar_flag, 23},
258 #ifdef ADJUST_GROUP_PRINT
259  {adjust_group_write_option.key, required_argument, &petar_flag, 24},
260 #endif
261  {nstep_dt_soft_kepler.key, required_argument, &petar_flag, 25},
262  {"help", no_argument, 0, 'h'},
263  {0,0,0,0}
264  };
265 
266  int opt_used = opt_used_pre;
267  int copt;
268  int option_index;
269  optind = 0; // reset getopt
270  while ((copt = getopt_long(argc, argv, "-i:a:t:s:o:r:b:n:G:u:T:f:p:w:h", long_options, &option_index)) != -1)
271  switch (copt) {
272  case 0:
273  switch (petar_flag) {
274 #ifdef ORBIT_SAMPLING
275  case 0:
276  n_split.value = atoi(optarg);
277  if(print_flag) n_split.print(std::cout);
278  opt_used += 2;
279  assert(n_split.value>=0);
280  break;
281 #endif
282  case 1:
283  search_vel_factor.value = atof(optarg);
284  if(print_flag) search_vel_factor.print(std::cout);
285  opt_used += 2;
286  update_rsearch_flag = true;
287  assert(search_vel_factor.value>0.0);
288  break;
289  case 2:
290  dt_limit_hard_factor.value = atof(optarg);
291  if(print_flag) dt_limit_hard_factor.print(std::cout);
292  opt_used += 2;
293  assert(dt_limit_hard_factor.value > 0.0);
294  break;
295  case 3:
296  dt_min_hermite_index.value = atoi(optarg);
297  if(print_flag) dt_min_hermite_index.print(std::cout);
298  opt_used += 2;
299  assert(dt_min_hermite_index.value > 0);
300  break;
301  case 4:
302  n_group_limit.value = atoi(optarg);
303  if(print_flag) n_group_limit.print(std::cout);
304  opt_used += 2;
305  assert(n_group_limit.value>0);
306  break;
307  case 5:
308  n_leaf_limit.value = atoi(optarg);
309  if(print_flag) n_leaf_limit.print(std::cout);
310  opt_used += 2;
311  assert(n_leaf_limit.value>0);
312  break;
313  case 6:
314  n_smp_ave.value = atoi(optarg);
315  if(print_flag) n_smp_ave.print(std::cout);
316  opt_used += 2;
317  assert(n_smp_ave.value>0.0);
318  break;
319  case 7:
320  e_err_ar.value = atof(optarg);
321  if(print_flag) e_err_ar.print(std::cout);
322  opt_used += 2;
323  assert(e_err_ar.value > 0.0);
324  break;
325  case 8:
326  eps.value = atof(optarg);
327  if(print_flag) eps.print(std::cout);
328  opt_used += 2;
329  assert(eps.value>=0.0);
330  break;
331  case 9:
332  sd_factor.value = atof(optarg);
333  if(print_flag) sd_factor.print(std::cout);
334  opt_used += 2;
335  assert(sd_factor.value>0.0);
336  break;
337  case 10:
338  ratio_r_cut.value = atof(optarg);
339  if(print_flag) ratio_r_cut.print(std::cout);
340  update_changeover_flag = true;
341  opt_used += 2;
342  assert(ratio_r_cut.value>0.0);
343  assert(ratio_r_cut.value<1.0);
344  break;
345  case 11:
346  r_bin.value = atof(optarg);
347  if(print_flag) r_bin.print(std::cout);
348  opt_used += 2;
349  assert(r_bin.value>=0.0);
350  break;
351  case 12:
352  search_peri_factor.value = atof(optarg);
353  if(print_flag) search_peri_factor.print(std::cout);
354  opt_used += 2;
355  assert(search_peri_factor.value>=1.0);
356  break;
357  case 13:
358  eta.value = atof(optarg);
359  if(print_flag) eta.print(std::cout);
360  opt_used += 2;
361  assert(eta.value>0.0);
362  break;
363 #ifdef HARD_CHECK_ENERGY
364  case 14:
365  e_err_hard.value = atof(optarg);
366  if(print_flag) e_err_hard.print(std::cout);
367  opt_used += 2;
368  break;
369 #endif
370  case 15:
371  step_limit_ar.value = atoi(optarg);
372  if(print_flag) step_limit_ar.print(std::cout);
373  opt_used += 2;
374  break;
375  case 16:
376  print_flag = false;
377  opt_used ++;
378  break;
379  case 17:
380  n_interrupt_limit.value = atoi(optarg);
381  if(print_flag) n_interrupt_limit.print(std::cout);
382  opt_used += 2;
383  assert(n_interrupt_limit.value>0);
384  break;
385  case 18:
386  interrupt_detection_option.value = atoi(optarg);
387  if(print_flag) interrupt_detection_option.print(std::cout);
388  opt_used += 2;
389  break;
390  case 19:
391  n_step_per_orbit.value = atof(optarg);
392  if(print_flag) n_step_per_orbit.print(std::cout);
393  opt_used += 2;
394  assert(n_step_per_orbit.value>=1.0);
395  break;
396 #ifdef STELLAR_EVOLUTION
397  case 20:
398  stellar_evolution_option.value = atoi(optarg);
399  if(print_flag) stellar_evolution_option.print(std::cout);
400  opt_used += 2;
401  break;
402 #endif
403  case 21:
404  r_escape.value = atof(optarg);
405  if(print_flag) r_escape.print(std::cout);
406  opt_used += 2;
407  break;
408  case 22:
409  r_search_min.value = atof(optarg);
410  if(print_flag) r_search_min.print(std::cout);
411  update_rsearch_flag = true;
412  opt_used += 2;
413  break;
414  case 23:
415  id_offset.value = atoi(optarg);
416  if(print_flag) id_offset.print(std::cout);
417  opt_used += 2;
418  break;
419 #ifdef ADJUST_GROUP_PRINT
420  case 24:
421  adjust_group_write_option.value = atoi(optarg);
422  if(print_flag) adjust_group_write_option.print(std::cout);
423  opt_used += 2;
424  break;
425 #endif
426  case 25:
427  nstep_dt_soft_kepler.value = atof(optarg);
428  if(print_flag) nstep_dt_soft_kepler.print(std::cout);
429  opt_used += 2;
430  break;
431  default:
432  break;
433  }
434  break;
435  case 'i':
436  data_format.value = atoi(optarg);
437  if(print_flag) data_format.print(std::cout);
438  assert(data_format.value>=0&&data_format.value<=3);
439  opt_used += 2;
440  break;
441  case 'a':
442  append_switcher.value=atoi(optarg);
443  if(print_flag) append_switcher.print(std::cout);
444  assert(append_switcher.value>=0&&append_switcher.value<=1);
445  opt_used += 2;
446  break;
447  case 't':
448  time_end.value = atof(optarg);
449  if(print_flag) time_end.print(std::cout);
450  opt_used += 2;
451  assert(time_end.value>=0.0);
452  break;
453  case 's':
454  dt_soft.value = atof(optarg);
455  if(print_flag) dt_soft.print(std::cout);
456  update_rsearch_flag = true;
457  opt_used += 2;
458  assert(dt_soft.value>=0.0);
459  break;
460  case 'o':
461  dt_snap.value = atof(optarg);
462  if(print_flag) dt_snap.print(std::cout);
463  opt_used += 2;
464  assert(dt_snap.value>0.0);
465  break;
466  case 'r':
467  r_out.value = atof(optarg);
468  if(print_flag) r_out.print(std::cout);
469  update_rsearch_flag = true;
470  update_changeover_flag = true;
471  opt_used += 2;
472  assert(r_out.value>=0.0);
473  break;
474  case 'b':
475  n_bin.value = atoi(optarg);
476  if(print_flag) n_bin.print(std::cout);
477  opt_used += 2;
478  assert(n_bin.value>=0);
479  break;
480  case 'n':
481  n_glb.value = atol(optarg);
482  if(print_flag) n_glb.print(std::cout);
483  opt_used += 2;
484  assert(n_glb.value>0);
485  break;
486  case 'G':
487  gravitational_constant.value = atof(optarg);
488  if(print_flag) gravitational_constant.print(std::cout);
489  opt_used += 2;
490  assert(gravitational_constant.value>0.0);
491  break;
492  case 'u':
493  unit_set.value = atoi(optarg);
494  if(print_flag) unit_set.print(std::cout);
495  opt_used += 2;
496  assert(unit_set.value>=0);
497  break;
498  case 'T':
499  theta.value = atof(optarg);
500  if(print_flag) theta.print(std::cout);
501  opt_used += 2;
502  assert(theta.value>=0.0);
503  break;
504  case 'f':
505  fname_snp.value = optarg;
506  if(print_flag) fname_snp.print(std::cout);
507  opt_used += 2;
508  break;
509  case 'p':
510  fname_par.value = optarg;
511  if(print_flag) {
512  fname_par.print(std::cout);
513  FILE* fpar_in;
514  if( (fpar_in = fopen(fname_par.value.c_str(),"r")) == NULL) {
515  fprintf(stderr,"Error: Cannot open file %s.\n", fname_par.value.c_str());
516  abort();
517  }
518  input_par_store.readAscii(fpar_in);
519  fclose(fpar_in);
520  }
521  opt_used += 2;
522 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
523  input_par_store.mpi_broadcast();
524  PS::Comm::barrier();
525 #endif
526  break;
527  case 'w':
528  write_style.value = atoi(optarg);
529  if(print_flag) write_style.print(std::cout);
530  opt_used += 2;
531  break;
532  case 'h':
533  if(print_flag){
534  std::cout<<"Usage: petar [option] filename"<<std::endl;
535  std::cout<<" filename: initial or restart data (snapshot) filename\n"
536  <<" Data file content:\n"
537  <<" First line: \n"
538  <<" 1. File_ID: 0 for initialization, else for restarting\n"
539  <<" 2. N_particle \n"
540  <<" 3. Time\n"
541 #ifdef RECORD_CM_IN_HEADER
542  <<" 4-6: offset of system centeral position (assuming galactic center is the coordiante origin in the case of Galpy)\n"
543  <<" 7-9: offset of system centeral velocity\n"
544 #endif
545  <<" Following lines:\n";
546  FPSoft::printTitleWithMeaning(std::cout,0,13);
547  std::cout<<" PS: (*) show initialization values which should be used together with FILE_ID = 0"<<std::endl;
548  std::cout<<" [formatted] indicates that the value is only for save, cannot be directly read"<<std::endl;
549  std::cout<<"Options:\n";
550  input_par_store.printHelp(std::cout, 2, 10, 23);
551  std::cout<<" --disable-print-info: "<<"Do not print information"<<std::endl;
552  std::cout<<" --disable-write-info: "<<"Do not write information"<<std::endl;
553  std::cout<<" -h(--help): print help"<<std::endl;
554  std::cout<<"*** PS: dt_soft: tree time step\n"
555  <<" r_in : transit function inner boundary radius\n"
556  <<" r_out: transit function outer boundary radius\n"
557  <<" sigma: half-mass radius velocity dispersion\n"
558  <<" n_bin: number of primordial binaries\n"
559  <<" <m> : averaged mass"<<std::endl;
560  }
561  return -1;
562  case '?':
563  opt_used +=2;
564  break;
565  default:
566  break;
567  }
568 
569  // count used options
570  opt_used ++;
571  //std::cout<<"Opt used:"<<opt_used<<std::endl;
572  if (opt_used<argc) {
573  fname_inp.value =argv[argc-1];
574  if(print_flag) std::cout<<"Reading data file name: "<<fname_inp.value<<std::endl;
575  }
576 
577  if(print_flag) std::cout<<"----- Finish reading input options of PeTar -----\n";
578 
579  return opt_used-1;
580  }
581 
583  bool checkParams() {
584 #ifdef ORBIT_SAMPLING
585  assert(n_split.value>=0);
586 #endif
587  assert(search_vel_factor.value>0.0);
588  assert(dt_limit_hard_factor.value > 0.0);
589  assert(dt_min_hermite_index.value > 0);
590  assert(e_err_ar.value > 0.0);
591  assert(eps.value>=0.0 && eps.value<=ratio_r_cut.value); // avoid incorrect self-potential correction after tree force, when eps>r_out, self-potential is G m /r_eps instead of G m/r_cut;
592  assert(sd_factor.value>0.0);
593  assert(ratio_r_cut.value>0.0);
594  assert(ratio_r_cut.value<1.0);
595  assert(r_bin.value>=0.0);
596  assert(search_peri_factor.value>=1.0);
597  assert(data_format.value>=0||data_format.value<=3);
598  assert(time_end.value>=0.0);
599  assert(dt_soft.value>=0.0);
600  assert(dt_snap.value>0.0);
601  assert(r_out.value>=0.0);
602  assert(n_bin.value>=0);
603  assert(n_glb.value>0);
604  assert(n_group_limit.value>0);
605  assert(n_interrupt_limit.value>0);
606  assert(n_leaf_limit.value>0);
607  assert(n_smp_ave.value>0.0);
608  assert(theta.value>=0.0);
609  assert(eta.value>0.0);
610  return true;
611  }
612 
613 };
614 
615 
617 class PeTar {
618 public:
619 #ifdef USE_QUAD
620  typedef PS::TreeForForceLong<ForceSoft, EPISoft, EPJSoft>::QuadrupoleWithSymmetrySearch TreeForce;
621 #else
622  typedef PS::TreeForForceLong<ForceSoft, EPISoft, EPJSoft>::MonopoleWithSymmetrySearch TreeForce;
623 #endif
624  typedef PS::ParticleSystem<FPSoft> SystemSoft;
625 
626  // For neighbor searching
627  typedef PS::TreeForForceShort<ForceSoft, EPISoft, EPJSoft>::Symmetry TreeNB;
628 
629  // IO
630  IOParamsPeTar input_parameters;
631 #ifdef BSE_BASE
632  IOParamsBSE bse_parameters;
633  std::string fbse_par_suffix = BSEManager::getBSEOutputFilenameSuffix();
634  std::string fsse_par_suffix = BSEManager::getSSEOutputFilenameSuffix();
635 #endif // BSE_BASE
636 #ifdef GALPY
637  IOParamsGalpy galpy_parameters;
638 #endif
639 
640 #ifdef PROFILE
641  PS::S32 dn_loop;
642  SysProfile profile;
643  SysCounts n_count;
644  SysCounts n_count_sum;
645  FDPSProfile tree_soft_profile;
646  FDPSProfile tree_nb_profile;
647  std::ofstream fprofile;
648 #endif
649 
650  Status stat;
651  std::ofstream fstatus;
652  PS::F64 time_kick;
653 
654  // escaper
655  Escaper escaper;
656  std::ofstream fesc;
657 
658  // file system
659  FileHeader file_header;
660  SystemSoft system_soft;
661 
662  // particle index map
663  std::map<PS::S64, PS::S32> id_adr_map;
664 
665  // domain
666  PS::S64 n_loop; // count for domain decomposition
667  PS::F64 domain_decompose_weight;
668  PS::DomainInfo dinfo;
669  PS::F64ort * pos_domain;
670 #ifdef FDPS_COMM
671  PS::CommInfo comm_info; // MPI communicator
672 #endif
673 
674  // tree time step manager
675  KickDriftStep dt_manager;
676 
677  // tree
678  TreeNB tree_nb;
679  TreeForce tree_soft;
680 
681 #ifdef GALPY
682  GalpyManager galpy_manager;
683 #endif
684 
685  // hard integrator
686  HardManager hard_manager;
687  SystemHard system_hard_one_cluster;
688  SystemHard system_hard_isolated;
689 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
690  SystemHard system_hard_connected;
691 #endif
692  int n_interrupt_glb;
693 
694  // mass change particle list
695  PS::ReallocatableArray<PS::S32> mass_modify_list;
696 
697  // remove list
698  PS::ReallocatableArray<PS::S32> remove_list;
699  PS::ReallocatableArray<PS::S64> remove_id_record;
700  PS::S32 n_remove;
701 
702  // search cluster
703  SearchCluster search_cluster;
704 
705  // safety check flag
706  bool read_parameters_flag;
707  bool read_data_flag;
708  bool initial_parameters_flag;
709  bool initial_step_flag;
710 
711  // MPI
712  PS::S32 my_rank;
713  PS::S32 n_proc;
714 
715  // FPDS flag
716  static bool initial_fdps_flag;
717 
719  PeTar():
720  input_parameters(),
721 #ifdef BSE_BASE
722  bse_parameters(),
723 #endif
724 #ifdef GALPY
725  galpy_parameters(),
726 #endif
727 #ifdef PROFILE
728  // profile
729  dn_loop(0), profile(), n_count(), n_count_sum(), tree_soft_profile(), fprofile(),
730 #endif
731  stat(), fstatus(), time_kick(0.0),
732  escaper(), fesc(),
733  file_header(), system_soft(), id_adr_map(),
734  n_loop(0), domain_decompose_weight(1.0), dinfo(), pos_domain(NULL),
735  dt_manager(),
736  tree_nb(), tree_soft(),
737 #ifdef GALPY
738  galpy_manager(),
739 #endif
740  hard_manager(), system_hard_one_cluster(), system_hard_isolated(),
741 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
742  system_hard_connected(),
743 #endif
744  n_interrupt_glb(0),
745  mass_modify_list(), remove_list(), remove_id_record(),
746  search_cluster(),
747  read_parameters_flag(false), read_data_flag(false), initial_parameters_flag(false), initial_step_flag(false) {
748  assert(initial_fdps_flag);
749  my_rank = PS::Comm::getRank();
750  n_proc = PS::Comm::getNumberOfProc();
751  }
752 
753 
755  PS::F64 regularTimeStep(const PS::F64 _dt) {
756  // regularize dt_tree
757  PS::F64 dt = 1.0;
758  if (_dt<1) while (dt>_dt) dt *= 0.5;
759  else {
760  while (dt<=_dt) dt *= 2.0;
761  dt *= 0.5;
762  }
763  return dt;
764  }
765 
767  void treeNeighborSearch() {
768 #ifdef PROFILE
769  profile.tree_nb.start();
770  tree_nb.clearNumberOfInteraction();
771  tree_nb.clearTimeProfile();
772 #endif
773 #ifdef USE_SIMD
774  tree_nb.calcForceAllAndWriteBack(SearchNeighborEpEpSimd(), system_soft, dinfo);
775 #elif USE_FUGAKU
776  tree_nb.calcForceAllAndWriteBack(SearchNeighborEpEpFugaku(), system_soft, dinfo);
777 #else
778  tree_nb.calcForceAllAndWriteBack(SearchNeighborEpEpNoSimd(), system_soft, dinfo);
779 #endif
780 
781 #ifdef PROFILE
782  tree_nb_profile += tree_nb.getTimeProfile();
783  //profile.tree_nb.barrier();
784  //PS::Comm::barrier();
785  profile.tree_nb.end();
786  profile.tree_nb.tbar = profile.tree_nb.time - tree_nb_profile.getTotalTime();
787 #endif
788  }
789 
791  void searchCluster() {
792 #ifdef PROFILE
793  profile.search_cluster.start();
794 #endif
795  // >2.1 search clusters ----------------------------------------
796  search_cluster.searchNeighborOMP<SystemSoft, TreeNB, EPJSoft>
797  (system_soft, tree_nb, pos_domain, 1.0, input_parameters.search_peri_factor.value);
798 
799  search_cluster.searchClusterLocal();
800  search_cluster.setIdClusterLocal();
801 
802  // >2.2 Send/receive connect cluster
803 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
804  search_cluster.connectNodes(pos_domain,tree_nb);
805  search_cluster.setIdClusterGlobalIteration();
806  search_cluster.sendAndRecvCluster(system_soft);
807 #endif
808 
809 #ifdef PROFILE
810  profile.search_cluster.barrier();
811  PS::Comm::barrier();
812  profile.search_cluster.end();
813 #endif
814  }
815 
817 
819  void createGroup(const PS::F64 _dt_tree) {
820 #ifdef PROFILE
821  profile.create_group.start();
822 #endif
823 
824  // >2.3 Find ARC groups and create artificial particles
825  // Set local ptcl_hard for isolated clusters
826  system_hard_isolated.setPtclForIsolatedMultiClusterOMP(system_soft, search_cluster.adr_sys_multi_cluster_isolated_, search_cluster.n_ptcl_in_multi_cluster_isolated_);
827 
828 //#ifdef CLUSTER_DEBUG
829 // for (PS::S32 i=0; i<n_loc; i++) system_soft[i].status = -1000000;
830 //#endif
831 
832  // Find groups and add artificial particles to global particle system
833  system_hard_isolated.findGroupsAndCreateArtificialParticlesOMP<SystemSoft, FPSoft>(system_soft, _dt_tree);
834 
835 //#ifdef CLUSTER_DEBUG
836 // not correct check, isolated clusters are not reset above
837 // for (PS::S32 i=0; i<n_loc; i++) {
838 // if (system_soft[i].status ==-1000000) {
839 // std::cerr<<"Find un initialized status i="<<i<<" status="<<system_soft[i].status<<std::endl;
840 // }
841 // }
842 //#endif
843  // update n_loc_iso, n_glb_iso for isolated clusters
844  //PS::S64 n_loc_iso = system_soft.getNumberOfParticleLocal();
845  //PS::S64 n_glb_iso = system_soft.getNumberOfParticleGlobal();
846 
847 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
848  // For connected clusters
849  system_hard_connected.setPtclForConnectedCluster(system_soft, search_cluster.mediator_sorted_id_cluster_, search_cluster.ptcl_recv_);
850  system_hard_connected.findGroupsAndCreateArtificialParticlesOMP<SystemSoft, FPSoft>(system_soft, _dt_tree);
851  // send updated particle back to original (set zero mass particle to origin)
852  search_cluster.writeAndSendBackPtcl(system_soft, system_hard_connected.getPtcl(), mass_modify_list);
853  mass_modify_list.resizeNoInitialize(0);
854  system_hard_connected.updateTimeWriteBack();
855 #endif
856 
857  // update total particle number including artificial particles
858  stat.n_all_loc = system_soft.getNumberOfParticleLocal();
859  stat.n_all_glb = system_soft.getNumberOfParticleGlobal();
860 
861  // >2.4 set adr/rank for artificial particles in GPS
862 #pragma omp parallel for
863  for(PS::S32 i=stat.n_real_loc; i<stat.n_all_loc; i++){
864  system_soft[i].rank_org = my_rank;
865  system_soft[i].adr = i;
866 #ifdef HARD_DEBUG
867  assert(system_soft[i].group_data.artificial.isArtificial());
868 #endif
869  }
870  // >3 Tree for force ----------------------------------------
871 #ifdef PROFILE
872  profile.create_group.barrier();
873  PS::Comm::barrier();
874  profile.create_group.end();
875 #endif
876  }
877 
879  void treeSoftForce() {
880 #ifdef PROFILE
881  profile.tree_soft.start();
882 
883  tree_soft.clearNumberOfInteraction();
884  tree_soft.clearTimeProfile();
885 #endif
886 
887 #ifdef USE_GPU
888  const PS::S32 n_walk_limit = 200;
889  const PS::S32 tag_max = 1;
890  PS::F64 eps2 = EPISoft::eps*EPISoft::eps;
891  PS::F64 rout2 = EPISoft::r_out*EPISoft::r_out;
892  PS::F64 G= ForceSoft::grav_const;
893 #ifdef PARTICLE_SIMULATOR_GPU_MULIT_WALK_INDEX
894  tree_soft.calcForceAllAndWriteBackMultiWalkIndex(CalcForceWithLinearCutoffCUDAMultiWalk(my_rank, eps2, rout2, G),
895  RetrieveForceCUDA,
896  tag_max,
897  system_soft,
898  dinfo,
899  n_walk_limit);
900 #else // no multi-walk index
901  tree_soft.calcForceAllAndWriteBackMultiWalk(CalcForceWithLinearCutoffCUDA(my_rank, eps2, rout2, G),
902  RetrieveForceCUDA,
903  tag_max,
904  system_soft,
905  dinfo,
906  n_walk_limit);
907 #endif // multi-walk index
908 
909 #elif USE_FUGAKU
910  PS::F64 eps2 = EPISoft::eps*EPISoft::eps;
911  PS::F64 rout2 = EPISoft::r_out*EPISoft::r_out;
912  PS::F64 G= ForceSoft::grav_const;
913  tree_soft.calcForceAllAndWriteBack(CalcForceEpEpWithLinearCutoffFugaku(eps2, rout2, G),
914 #ifdef USE_QUAD
915  CalcForceEpSpQuadFugaku(eps2, G),
916 #else // no quad
917  CalcForceEpSpMonoFugaku(eps2, G),
918 #endif // end quad
919  system_soft,
920  dinfo);
921 
922 #elif USE_SIMD // end use_gpu
923  tree_soft.calcForceAllAndWriteBack(CalcForceEpEpWithLinearCutoffSimd(),
924 #ifdef USE_QUAD
925  CalcForceEpSpQuadSimd(),
926 #else // no quad
927  CalcForceEpSpMonoSimd(),
928 #endif // end quad
929  system_soft,
930  dinfo);
931 #else // end use_simd
932  tree_soft.calcForceAllAndWriteBack(CalcForceEpEpWithLinearCutoffNoSimd(),
933 #ifdef USE_QUAD
934  CalcForceEpSpQuadNoSimd(),
935 #else
936  CalcForceEpSpMonoNoSimd(),
937 #endif
938  system_soft,
939  dinfo);
940 #endif // end else
941 
942 #ifdef PROFILE
943  n_count.ep_ep_interact += tree_soft.getNumberOfInteractionEPEPLocal();
944  n_count_sum.ep_ep_interact += tree_soft.getNumberOfInteractionEPEPGlobal();
945  n_count.ep_sp_interact += tree_soft.getNumberOfInteractionEPSPLocal();
946  n_count_sum.ep_sp_interact += tree_soft.getNumberOfInteractionEPSPGlobal();
947 
948  tree_soft_profile += tree_soft.getTimeProfile();
949  domain_decompose_weight = tree_soft_profile.calc_force;
950 
951  //profile.tree_soft.barrier();
952  //PS::Comm::barrier();
953  profile.tree_soft.end();
954  profile.tree_soft.tbar = profile.tree_soft.time - tree_soft_profile.getTotalTime();
955 #endif
956  }
957 
959  void treeForceCorrectChangeoverTreeNeighbor() {
960  // all particles
961  SystemHard::correctForceWithCutoffTreeNeighborOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, system_soft.getNumberOfParticleLocal(), hard_manager.ap_manager);
962  }
963 
965  void treeForceCorrectChangeover() {
966 #ifdef PROFILE
967  profile.force_correct.start();
968 #endif
969 
970 #ifdef CORRECT_FORCE_DEBUG
971  // backup particle data
972  PS::S64 n_loc_all = system_soft.getNumberOfParticleLocal();
973 
974  FPSoft psys_bk[n_loc_all];
975 #pragma omp parallel for
976  for (int i=0; i<n_loc_all; i++)
977  psys_bk[i] = system_soft[i];
978 #endif
979 
980  // single
981  system_hard_one_cluster.correctPotWithCutoffOMP(system_soft, search_cluster.getAdrSysOneCluster());
982 
983  // Isolated clusters
984  system_hard_isolated.correctForceWithCutoffClusterOMP(system_soft);
985 
986 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
987  // Connected clusters
988  system_hard_connected.correctForceWithCutoffTreeNeighborAndClusterOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, search_cluster.getAdrSysConnectClusterSend());
989 #endif
990 
991 #ifdef CORRECT_FORCE_DEBUG
992  // switch backup and sys particle data
993 #pragma omp parallel for
994  for (int i=0; i<n_loc_all; i++) {
995  FPSoft ptmp = psys_bk[i];
996  psys_bk[i] = system_soft[i];
997  system_soft[i] = ptmp;
998  }
999 
1000  // all particles
1001  SystemHard::correctForceWithCutoffTreeNeighborOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, n_loc, hard_manager.ap_manager);
1002 
1003  // single
1004  //system_hard_one_cluster.correctPotWithCutoffOMP(system_soft, search_cluster.getAdrSysOneCluster());
1005  // Isolated clusters
1006  //system_hard_isolated.correctForceWithCutoffTreeNeighborOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, n_loc);
1007 //#ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1008 // // Connected clusters
1009 // system_hard_connected.correctForceWithCutoffTreeNeighborOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, n_loc_all);
1010 //#endif
1011 
1012  for (int i=0; i<n_loc_all; i++) {
1013  PS::F64vec dacci=psys_bk[i].acc- system_soft[i].acc;
1014  PS::F64 dpoti = psys_bk[i].pot_tot- system_soft[i].pot_tot;
1015  if(dacci*dacci/(system_soft[i].acc*system_soft[i].acc)>1e-8) {
1016  std::cerr<<"Corrected Acc diff >1e-8: i "<<i<<" acc(tree): "<<system_soft[i].acc<<" acc(cluster): "<<psys_bk[i].acc<<std::endl;
1017  abort();
1018  }
1019  // Notice, the cluster members can include particles are not in the tree neighbor searching. If the extra neighbors in clusters are group members. The correct from clusters will be more accurate than tree neighbor correction. Since the potential is calculated from the real members instead of artificial particles. This can result in different potential
1020  if(abs(dpoti/system_soft[i].pot_tot)>1e-6) {
1021  std::cerr<<"Corrected pot diff >1e-6: i "<<i<<" pot(tree): "<<system_soft[i].pot_tot<<" pot(cluster): "<<psys_bk[i].pot_tot<<std::endl;
1022  abort();
1023  }
1024  }
1025 
1026 #endif
1027 #ifdef PROFILE
1028  profile.force_correct.barrier();
1029  PS::Comm::barrier();
1030  profile.force_correct.end();
1031 #endif
1032  }
1033 
1035  void externalForce() {
1036 #ifdef PROFILE
1037  profile.other.start();
1038 #endif
1039 
1040 #ifdef GALPY
1041  // external force and potential
1042  // update the types and arguments
1043  galpy_manager.updatePotential(stat.time, true);
1044 
1045  galpy_manager.resetPotAcc();
1046  galpy_manager.calcMovePotAccFromPot(stat.time, &stat.pcm.pos[0]);
1047 
1048  PS::S64 n_loc_all = system_soft.getNumberOfParticleLocal();
1049 #pragma omp parallel for
1050  for (int i=0; i<n_loc_all; i++) {
1051  auto& pi = system_soft[i];
1052  double acc[3], pot;
1053 #ifdef RECORD_CM_IN_HEADER
1054  PS::F64vec pos_correct=pi.pos + stat.pcm.pos;
1055  galpy_manager.calcAccPot(acc, pot, stat.time, input_parameters.gravitational_constant.value*pi.mass, &pos_correct[0], &pi.pos[0]);
1056 #else
1057  PS::F64vec pos_center=pi.pos - stat.pcm.pos;
1058  galpy_manager.calcAccPot(acc, pot, stat.time, input_parameters.gravitational_constant.value*pi.mass, &pi.pos[0], &pos_center[0]);
1059 #endif
1060  assert(!std::isinf(acc[0]));
1061  assert(!std::isnan(acc[0]));
1062  assert(!std::isinf(pot));
1063  assert(!std::isnan(pot));
1064  pi.acc[0] += acc[0];
1065  pi.acc[1] += acc[1];
1066  pi.acc[2] += acc[2];
1067  pi.pot_tot += pot;
1068  pi.pot_soft += pot;
1069 #ifdef EXTERNAL_POT_IN_PTCL
1070  pi.pot_ext = pot;
1071 #endif
1072  }
1073 #endif //GALPY
1074 
1075 #ifdef PROFILE
1076  profile.other.barrier();
1077  PS::Comm::barrier();
1078  profile.other.end();
1079 #endif
1080  }
1081 
1082  // correct c.m. vel due to external potential to avoid large rsearch
1083  void correctPtclVelCM(const PS::F64& _dt) {
1084  PS::F64vec dv=PS::F64vec(0.0);
1085 
1086 #ifdef GALPY
1087  PS::F64vec acc;
1088  PS::F64 pot;
1089  // evaluate center of mass acceleration
1090  PS::F64vec pos_zero=PS::F64vec(0.0);
1091  // set zero mass to avoid duplicate anti force to potential set
1092  galpy_manager.calcAccPot(&acc[0], pot, stat.time, 0, &stat.pcm.pos[0], &pos_zero[0]);
1093  dv = acc*_dt;
1094 #endif
1095 
1096  stat.pcm.vel += dv;
1097  for (int i=0; i<stat.n_all_loc; i++) system_soft[i].vel -= dv;
1098 
1099  //auto& adr = search_cluster.getAdrSysOneCluster();
1100  //for (int i=0; i<adr.size(); i++) system_soft[adr[i]].vel -= dv;
1101 
1102  auto& ptcl_iso=system_hard_isolated.getPtcl();
1103  for (int i=0; i<ptcl_iso.size(); i++) ptcl_iso[i].vel -= dv;
1104 
1105  //const PS::S64 n_tot= system_soft.getNumberOfParticleLocal();
1106  //const PS::S32 n_artifical_per_group = hard_manager.ap_manager.getArtificialParticleN();
1107  //for(PS::S32 i=stat.n_real_loc; i<n_tot; i+= n_artifical_per_group) {
1108  // auto* pcm = hard_manager.ap_manager.getCMParticles(&(system_soft[i]));
1109  // pcm->vel -= dv;
1110  //}
1111 
1112 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1113  auto& ptcl_con=system_hard_connected.getPtcl();
1114  for (int i=0; i<ptcl_con.size(); i++) ptcl_con[i].vel -= dv;
1115 #endif
1116 
1117  //Ptcl::vel_cm += dv;
1118 //#ifdef PETAR_DEBUG
1119 // if (my_rank==0) std::cerr<<"Ptcl::vel_cm: "<<Ptcl::vel_cm<<std::endl;
1120 //#endif
1121  }
1122 
1123 #ifdef KDKDK_4TH
1124  void GradientKick() {
1126 #ifdef PROFILE
1127  profile.tree_soft.start();
1128  tree_soft.clearNumberOfInteraction();
1129  tree_soft.clearTimeProfile();
1130 #endif
1131  // correction calculation
1132  //tree_soft.setParticaleLocalTree(system_soft, false);
1133 
1134  tree_soft.calcForceAllAndWriteBack(CalcCorrectEpEpWithLinearCutoffNoSimd(),
1135 #ifdef USE_QUAD
1136  CalcForceEpSpQuadNoSimd(),
1137 #else
1138  CalcForceEpSpMonoNoSimd(),
1139 #endif
1140  system_soft,
1141  dinfo);
1142 
1143 #ifdef PROFILE
1144  n_count.ep_ep_interact += tree_soft.getNumberOfInteractionEPEPLocal();
1145  n_count_sum.ep_ep_interact += tree_soft.getNumberOfInteractionEPEPGlobal();
1146  n_count.ep_sp_interact += tree_soft.getNumberOfInteractionEPSPLocal();
1147  n_count_sum.ep_sp_interact += tree_soft.getNumberOfInteractionEPSPGlobal();
1148 
1149  tree_soft_profile += tree_soft.getTimeProfile();
1150  domain_decompose_weight += tree_soft_profile.calc_force;
1151 
1152  profile.tree_soft.barrier();
1153  PS::Comm::barrier();
1154  profile.tree_soft.end();
1155  profile.force_correct.start();
1156 #endif
1157 
1158  // Isolated clusters
1159  system_hard_isolated.correctForceWithCutoffClusterOMP(system_soft, true);
1160 
1161 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1162  // Connected clusters
1163  system_hard_connected.correctForceWithCutoffTreeNeighborAndClusterOMP<SystemSoft, FPSoft, TreeForce, EPJSoft>(system_soft, tree_soft, search_cluster.getAdrSysConnectClusterSend(), true);
1164 #endif
1165 
1166 #ifdef PROFILE
1167  profile.force_correct.barrier();
1168  PS::Comm::barrier();
1169  profile.force_correct.end();
1170 #endif
1171 // if (true) {
1172 
1173 
1174 // for debug
1175 // FILE* fdump;
1176 // if ( (fdump = fopen("acorr.dump","w")) == NULL) {
1177 // fprintf(stderr,"Error: Cannot open file acorr.dump\n");
1178 // abort();
1179 // }
1180 // for (int i=0; i<n_loc; i++) {
1181 // fprintf(fdump, "%d %.12g %.12g %.12g %.12g %.12g %.12g\n", i, system_soft[i].acc[0], system_soft[i].acc[1], system_soft[i].acc[2], system_soft[i].acorr[0], system_soft[i].acorr[1], system_soft[i].acorr[2]);
1182 // }
1183 // fclose(fdump);
1184 // abort();
1185 // debug
1186 // }
1187 
1188  }
1189 #endif
1190 
1192  /* modify the velocity of particle in global system
1193  reset particle type to single
1194  @param[in,out] _sys: particle system
1195  @param[in]: _dt: tree step
1196  @param[in]; _adr: address for single particles
1197  */
1198  void kickOne(SystemSoft & _sys,
1199  const PS::F64 _dt,
1200  const PS::ReallocatableArray<PS::S32>& _adr) {
1201  const PS::S64 n= _adr.size();
1202 #pragma omp parallel for
1203  for(PS::S32 i=0; i<n; i++){
1204  const PS::S32 k=_adr[i];
1205 #ifdef KDKDK_4TH
1206  _sys[k].vel += _dt*(_sys[k].acc + 9.0/192.0*_dt*_dt*_sys[k].acorr);
1207 #else
1208  _sys[k].vel += _sys[k].acc * _dt;
1209 #endif
1210  _sys[k].group_data.artificial.setParticleTypeToSingle();
1211  }
1212  }
1213 
1215 
1221  template<class Tptcl>
1222  void kickClusterAndRecoverGroupMemberMass(SystemSoft& _sys,
1223  PS::ReallocatableArray<Tptcl>& _ptcl,
1224  const PS::F64 _dt) {
1225  assert(Ptcl::group_data_mode == GroupDataMode::artificial);
1226  const PS::S64 n= _ptcl.size();
1227 #pragma omp parallel for
1228  for(PS::S32 i=0; i<n; i++) {
1229  auto& pi_artificial = _ptcl[i].group_data.artificial;
1230  // if is group member, recover mass and kick due to c.m. force
1231  if (pi_artificial.isMember()) {
1232  const PS::S64 cm_adr = _ptcl[i].getParticleCMAddress();
1233  if (cm_adr>0) {
1234 #ifdef HARD_DEBUG
1235  assert(pi_artificial.getMassBackup()>0);
1236  assert(cm_adr<stat.n_all_loc);
1237 #endif
1238  _ptcl[i].mass = pi_artificial.getMassBackup();
1239 #ifdef KDKDK_4TH
1240  _ptcl[i].vel += _dt*(_sys[cm_adr].acc + 9.0/192.0*_dt*_dt*_sys[cm_adr].acorr);
1241 #else
1242 #ifdef NAN_CHECK_DEBUG
1243  assert(!std::isnan(_sys[cm_adr].acc[0]));
1244  assert(!std::isnan(_sys[cm_adr].acc[1]));
1245  assert(!std::isnan(_sys[cm_adr].acc[2]));
1246 #endif
1247  _ptcl[i].vel += _sys[cm_adr].acc * _dt;
1248 #endif
1249  continue;
1250  }
1251  }
1252  // non-member particle
1253  const PS::S64 i_adr =_ptcl[i].adr_org;
1254  if(i_adr>=0) {
1255  // not remote particles
1256 #ifdef HARD_DEBUG
1257  assert(pi_artificial.getStatus()==_sys[i_adr].group_data.artificial.getStatus());
1258 #endif
1259 
1260 #ifdef KDKDK_4TH
1261  _ptcl[i].vel += _dt*(_sys[i_adr].acc + 9.0/192.0*_dt*_dt*_sys[i_adr].acorr);
1262 #else
1263  _ptcl[i].vel += _sys[i_adr].acc * _dt;
1264 #endif
1265  }
1266  }
1267  }
1268 
1270 
1275  void kickSend(SystemSoft& _sys,
1276  const PS::ReallocatableArray<PS::S32>& _adr_ptcl_send,
1277  const PS::F64 _dt) {
1278  assert(Ptcl::group_data_mode == GroupDataMode::artificial);
1279  const PS::S64 n= _adr_ptcl_send.size();
1280 #pragma omp parallel for
1281  for(PS::S32 i=0; i<n; i++) {
1282  const PS::S64 adr = _adr_ptcl_send[i];
1283  // if it is group member with artificial particles, should not do kick since the required c.m. forces is on remote nodes;
1284  // if it is group member without artificial particles, also kick
1285  if(_sys[adr].group_data.artificial.isSingle() || (_sys[adr].group_data.artificial.isMember() && _sys[adr].getParticleCMAddress()<0)) {
1286  _sys[adr].vel += _sys[adr].acc * _dt;
1287 #ifdef KDKDK_4TH
1288  _sys[adr].vel += _dt*_dt* _sys[adr].acorr /48;
1289 #endif
1290  }
1291 
1292 #ifdef HARD_DEBUG
1293  if(_sys[adr].group_data.artificial.isSingle()) assert(_sys[adr].mass>0);
1294  else assert(_sys[adr].group_data.artificial.getMassBackup()>0);
1295 #endif
1296  }
1297  }
1298 
1300 
1306  void kickCM(SystemSoft& _sys,
1307  const PS::S32 _adr_artificial_start,
1308  ArtificialParticleManager& _ap_manager,
1309  const PS::F64 _dt) {
1310  const PS::S64 n_tot= _sys.getNumberOfParticleLocal();
1311  const PS::S32 n_artifical_per_group = _ap_manager.getArtificialParticleN();
1312 #pragma omp parallel for
1313  for(PS::S32 i=_adr_artificial_start; i<n_tot; i+= n_artifical_per_group) {
1314  auto* pcm = _ap_manager.getCMParticles(&(_sys[i]));
1315  pcm->vel += pcm->acc * _dt;
1316 #ifdef KDKDK_4TH
1317  pcm->vel += _dt*_dt* pcm->acorr /48;
1318 #endif
1319 #ifdef HARD_DEBUG
1320  assert(pcm->group_data.artificial.isCM());
1321 #endif
1322  }
1323  }
1324 
1326  void driftSingle(SystemSoft & system,
1327  const PS::F64 dt){
1328  const PS::S32 n = system.getNumberOfParticleLocal();
1329 #pragma omp parallel for
1330  for(PS::S32 i=0; i<n; i++){
1331  if(system[i].n_ngb <= 0){
1332  system[i].pos += system[i].vel * dt;
1333 //#ifdef STELLAR_EVOLUTION
1334  // shift time interrupt in order to get consistent time for stellar evolution in the next drift
1335  //system[i].time_record -= dt;
1336  //system[i].time_interrupt -= dt;
1337 //#endif
1338  }
1339  }
1340  }
1341 
1343  void kick(const PS::F64 _dt_kick) {
1344 #ifdef PROFILE
1345  profile.kick.start();
1346 #endif
1347 
1349  // single and reset particle type to single (due to binary disruption)
1350  kickOne(system_soft, _dt_kick, search_cluster.getAdrSysOneCluster());
1351  // isolated
1352  kickClusterAndRecoverGroupMemberMass(system_soft, system_hard_isolated.getPtcl(), _dt_kick);
1353  // c.m. artificial
1354  kickCM(system_soft, stat.n_real_loc, hard_manager.ap_manager, _dt_kick);
1355 
1356 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1357  // connected
1358  kickClusterAndRecoverGroupMemberMass(system_soft, system_hard_connected.getPtcl(), _dt_kick);
1359  // sending list for connected clusters
1360  kickSend(system_soft, search_cluster.getAdrSysConnectClusterSend(), _dt_kick);
1361  // send kicked particle from sending list, and receive remote single particle
1362  search_cluster.SendSinglePtcl(system_soft, system_hard_connected.getPtcl());
1363 #endif
1364 
1365 #ifdef GALPY
1366  galpy_manager.kickMovePot(_dt_kick);
1367 #endif
1368 
1369 #ifdef RECORD_CM_IN_HEADER
1370  // correct Ptcl:vel_cm
1371  correctPtclVelCM(_dt_kick);
1372 #endif
1373 
1374 
1375  time_kick += _dt_kick;
1376 
1377 #ifdef HARD_DEBUG
1378  PS::S32 kick_regist[stat.n_real_loc];
1379  for(int i=0; i<stat.n_real_loc; i++) kick_regist[i] = 0;
1380  for(int i=0; i<search_cluster.getAdrSysOneCluster().size(); i++) {
1381  kick_regist[search_cluster.getAdrSysOneCluster()[i]]++;
1382  }
1383  for(int i=0; i<system_hard_isolated.getPtcl().size(); i++) {
1384  PS::S64 adr= system_hard_isolated.getPtcl()[i].adr_org;
1385  assert(adr>=0);
1386  kick_regist[adr]++;
1387  }
1388 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1389  for(int i=0; i<system_hard_connected.getPtcl().size(); i++) {
1390  PS::S64 adr= system_hard_connected.getPtcl()[i].adr_org;
1391  if(adr>=0) kick_regist[adr]++;
1392  }
1393  for(int i=0; i<search_cluster.getAdrSysConnectClusterSend().size(); i++) {
1394  PS::S64 adr= search_cluster.getAdrSysConnectClusterSend()[i];
1395  kick_regist[adr]++;
1396  }
1397 #endif
1398  for(int i=0; i<stat.n_real_loc; i++) {
1399  assert(kick_regist[i]==1);
1400  }
1401 #endif
1402 
1403 #ifdef PROFILE
1404  profile.kick.barrier();
1405  PS::Comm::barrier();
1406  profile.kick.end();
1407 #endif
1408  }
1409 
1411 
1413  PS::S32 drift(const PS::F64 _dt_drift) {
1415  //system_hard_one_cluster.setTimeOrigin(stat.time);
1416  //system_hard_isolated.setTimeOrigin(stat.time);
1417  //system_hard_connected.setTimeOrigin(stat.time);
1419 
1420 #ifdef PROFILE
1421  profile.hard_single.start();
1422 #endif
1423  system_hard_one_cluster.initializeForOneCluster(search_cluster.getAdrSysOneCluster().size());
1425  system_hard_one_cluster.setPtclForOneClusterOMP(system_soft, search_cluster.getAdrSysOneCluster());
1426  system_hard_one_cluster.driveForOneClusterOMP(_dt_drift);
1427  //system_hard_one_cluster.writeBackPtclForOneClusterOMP(system_soft, search_cluster.getAdrSysOneCluster());
1428  system_hard_one_cluster.writeBackPtclForOneClusterOMP(system_soft, mass_modify_list);
1430 
1431 #ifdef PROFILE
1432  profile.hard_single.barrier();
1433  PS::Comm::barrier();
1434  profile.hard_single.end();
1435 #endif
1436 
1439 #ifdef PROFILE
1440  profile.hard_isolated.start();
1441 #endif
1442  // reset slowdown energy correction
1443  system_hard_isolated.energy.resetEnergyCorrection();
1444  // integrate multi cluster A
1445  system_hard_isolated.driveForMultiClusterOMP(_dt_drift, &(system_soft[0]));
1446  //system_hard_isolated.writeBackPtclForMultiCluster(system_soft, search_cluster.adr_sys_multi_cluster_isolated_,remove_list);
1447  PS::S32 n_interrupt_isolated = system_hard_isolated.getNumberOfInterruptClusters();
1448  if(n_interrupt_isolated==0) system_hard_isolated.writeBackPtclForMultiCluster(system_soft, mass_modify_list);
1449  // integrate multi cluster A
1450 
1451 #ifdef PROFILE
1452  n_count.hard_interrupt += n_interrupt_isolated;
1453  profile.hard_isolated.barrier();
1454 #endif
1455 
1456 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1457  n_interrupt_glb = PS::Comm::getSum(n_interrupt_isolated);
1458 #else
1459  n_interrupt_glb = n_interrupt_isolated;
1460 #endif
1461 
1462 #ifdef PROFILE
1463  n_count_sum.hard_interrupt += n_interrupt_glb;
1464  profile.hard_isolated.end();
1465 #endif
1466 
1468 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1469 #ifdef PROFILE
1471  profile.hard_connected.start();
1472 #endif
1473  // reset slowdown energy correction
1474  system_hard_connected.energy.resetEnergyCorrection();
1475  // integrate multi cluster B
1476  system_hard_connected.driveForMultiClusterOMP(_dt_drift, &(system_soft[0]));
1477  PS::S32 n_interrupt_connected = system_hard_connected.getNumberOfInterruptClusters();
1478 
1479 #ifdef PROFILE
1480  n_count.hard_interrupt += n_interrupt_connected;
1481  profile.hard_connected.barrier();
1482 #endif
1483 
1484  PS::S32 n_interrupt_connected_glb = PS::Comm::getSum(n_interrupt_connected);
1485 
1486 #ifdef PROFILE
1487  n_count_sum.hard_interrupt += n_interrupt_connected_glb;
1488  profile.hard_connected.end();
1489  profile.hard_connected.start();
1490 #endif
1491 
1492 
1493  if (n_interrupt_connected_glb==0) {
1494  search_cluster.writeAndSendBackPtcl(system_soft, system_hard_connected.getPtcl(), mass_modify_list);
1495  system_hard_connected.updateTimeWriteBack();
1496  }
1497  // integrate multi cluster B
1498 
1499  n_interrupt_glb += n_interrupt_connected_glb;
1500 #ifdef PROFILE
1501  profile.hard_connected.barrier();
1502  PS::Comm::barrier();
1503  profile.hard_connected.end();
1504 #endif
1505 #endif
1506  // drift cm
1507  stat.pcm.pos += stat.pcm.vel*_dt_drift;
1508 
1509 #ifdef GALPY
1510  galpy_manager.driftMovePot(_dt_drift);
1511 #endif
1512 
1513  if (n_interrupt_glb==0) Ptcl::group_data_mode = GroupDataMode::cm;
1514 
1515  return n_interrupt_glb;
1516  }
1517 
1519 
1521  PS::S32 finishInterruptDrift() {
1522 #ifdef PROFILE
1523  profile.hard_interrupt.start();
1524  profile.hard_isolated.start();
1525 #endif
1526 
1527  // finish interrupt clusters first
1528  // isolated clusters
1529  PS::S32 n_interrupt_isolated = 0;
1530  if (system_hard_isolated.getNumberOfInterruptClusters()>0) {
1531  system_hard_isolated.finishIntegrateInterruptClustersOMP();
1532  n_interrupt_isolated = system_hard_isolated.getNumberOfInterruptClusters();
1533  if(n_interrupt_isolated==0) system_hard_isolated.writeBackPtclForMultiCluster(system_soft, mass_modify_list);
1534  }
1535 
1536 #ifdef PROFILE
1537  n_count.hard_interrupt += n_interrupt_isolated;
1538  profile.hard_isolated.barrier();
1539 #endif
1540 
1541 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1542  n_interrupt_glb = PS::Comm::getSum(n_interrupt_isolated);
1543 #else
1544  n_interrupt_glb = n_interrupt_isolated;
1545 #endif
1546 
1547 #ifdef PROFILE
1548  n_count_sum.hard_interrupt += n_interrupt_glb;
1549  profile.hard_isolated.end();
1550 #endif
1551 
1552 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1553 #ifdef PROFILE
1554  profile.hard_connected.start();
1555 #endif
1556 
1557  // connected clusters
1558  PS::S32 n_interrupt_connected = 0;
1559  if (system_hard_connected.getNumberOfInterruptClusters()>0) {
1560  system_hard_connected.finishIntegrateInterruptClustersOMP();
1561  n_interrupt_connected = system_hard_connected.getNumberOfInterruptClusters();
1562  }
1563 
1564 #ifdef PROFILE
1565  n_count.hard_interrupt += n_interrupt_connected;
1566  profile.hard_connected.barrier();
1567 #endif
1568 
1569  PS::S32 n_interrupt_connected_glb = PS::Comm::getSum(n_interrupt_connected);
1570 
1571 #ifdef PROFILE
1572  n_count_sum.hard_interrupt += n_interrupt_connected_glb;
1573  profile.hard_connected.end();
1574  profile.hard_connected.start();
1575 #endif
1576 
1577  if (n_interrupt_connected_glb==0) {
1578  PS::ReallocatableArray<PS::S32> mass_modify_list;
1579  search_cluster.writeAndSendBackPtcl(system_soft, system_hard_connected.getPtcl(), mass_modify_list);
1580  system_hard_connected.updateTimeWriteBack();
1581  }
1582 
1583  n_interrupt_glb += n_interrupt_connected_glb;
1584 
1585 #ifdef PROFILE
1586  profile.hard_connected.barrier();
1587  PS::Comm::barrier();
1588  profile.hard_connected.end();
1589 #endif
1590 #endif
1591 
1592 #ifdef HARD_INTERRUPT_PRINT
1593  if (n_interrupt_glb>0) {
1594  std::cerr<<"Interrupt detected, number: "<<n_interrupt_glb<<std::endl;
1595  }
1596 #endif
1597 
1598 #ifdef PROFILE
1599  profile.hard_interrupt.barrier();
1600  PS::Comm::barrier();
1601  profile.hard_interrupt.end();
1602 #endif
1603 
1604  if (n_interrupt_glb==0) Ptcl::group_data_mode = GroupDataMode::cm;
1605 
1606  // if interrupt cluster still exist, return the number immediately without new integration.
1607  return n_interrupt_glb;
1608  }
1609 
1610 
1612 
1618  inline PS::F64 calcDtLimit(const PS::F64 _time,
1619  const PS::F64 _dt_max,
1620  const PS::F64 _dt_min){
1621  // for first step, the maximum time step is OK
1622  if(_time==0.0) return _dt_max;
1623  else {
1624  // the binary tree for current time position in block step
1625  PS::U64 bitmap = _time/_dt_min;
1626 //#ifdef __GNUC__
1627 // PS::S64 dts = __builtin_ctz(bitmap) ;
1628 // PS::U64 c = (1<<dts);
1630 //#else
1631 
1632  // block step multiply factor
1633  PS::U64 c=1;
1634  // find the last zero in the binary tree to obtain the current block step level
1635  while((bitmap&1)==0) {
1636  bitmap = (bitmap>>1);
1637  c = (c<<1);
1638  }
1639 //#endif
1640  // return the maximum step allown
1641  return std::min(c*_dt_min,_dt_max);
1642  }
1643  }
1644 
1646  bool checkTimeConsistence() {
1647  if (abs(time_kick-stat.time)>1e-13) {
1648  std::cerr<<"Error: kick time ("<<time_kick<<") and system time ("<<stat.time<<") are inconsistent!"<<std::endl;
1649  abort();
1650  }
1651  time_kick = stat.time; // escape the problem of round-off error
1652  if (stat.n_real_glb>1) {
1653  assert(stat.time == system_hard_one_cluster.getTimeOrigin());
1654  assert(stat.time == system_hard_isolated.getTimeOrigin());
1655 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1656  assert(stat.time == system_hard_connected.getTimeOrigin());
1657 #endif
1658  }
1659  return true;
1660  }
1661 
1663  void writeBackHardParticles() {
1664 #ifdef PROFILE
1665  profile.output.start();
1666 #endif
1667  system_hard_isolated.writeBackPtclForMultiCluster(system_soft, mass_modify_list);
1668 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1669  // update gloabl particle system and send receive remote particles
1670  search_cluster.writeAndSendBackPtcl(system_soft, system_hard_connected.getPtcl(), mass_modify_list);
1671  system_hard_connected.updateTimeWriteBack();
1672 #endif
1673  mass_modify_list.resizeNoInitialize(0);
1674 
1675 #ifdef PROFILE
1676  profile.output.barrier();
1677  PS::Comm::barrier();
1678  profile.output.end();
1679 #endif
1680  }
1681 
1682  // update group_data.cm to pcm data for search cluster after restart
1683  void setParticleGroupDataToCMData() {
1684 #ifdef PROFILE
1685  profile.search_cluster.start();
1686 #endif
1687 #ifdef CLUSTER_VELOCITY
1688  // update status and mass_bk to pcm data for search cluster after restart
1689  system_hard_one_cluster.resetParticleGroupData(system_soft);
1690  system_hard_isolated.setParticleGroupDataToCMData(system_soft);
1691 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1692  system_hard_connected.setParticleGroupDataToCMData(system_soft);
1693  search_cluster.writeAndSendBackPtcl(system_soft, system_hard_connected.getPtcl(), mass_modify_list);
1694  system_hard_connected.updateTimeWriteBack();
1695  mass_modify_list.resizeNoInitialize(0);
1696 #endif
1697 #ifdef PROFILE
1698  profile.search_cluster.barrier();
1699  PS::Comm::barrier();
1700  profile.search_cluster.end();
1701 #endif
1702  }
1703 
1704 #endif
1705 
1707  void correctForceChangeOverUpdate() {
1708 #ifdef PROFILE
1709  profile.force_correct.start();
1710 #endif
1711  // correct changeover for first step
1712  // Isolated clusters
1713  system_hard_isolated.correctForceForChangeOverUpdateOMP<SystemSoft, TreeForce, EPJSoft>(system_soft, tree_soft);
1714 
1715 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1716  // Connected clusters
1717  auto& adr_send = search_cluster.getAdrSysConnectClusterSend();
1718  system_hard_connected.correctForceForChangeOverUpdateOMP<SystemSoft, TreeForce, EPJSoft>(system_soft, tree_soft, adr_send.getPointer(), adr_send.size());
1719 #endif
1720 
1721 #ifdef PROFILE
1722  profile.force_correct.barrier();
1723  PS::Comm::barrier();
1724  profile.force_correct.end();
1725 #endif
1726 
1727  }
1728 
1730 
1733  void domainDecompose(const bool _enforce=false) {
1734 #ifdef PROFILE
1735  // > 6. Domain decomposition
1736  profile.domain.start();
1737 #endif
1738  // Domain decomposition, parrticle exchange and force calculation
1739  if(n_loop % 16 == 0 || _enforce) {
1740  dinfo.decomposeDomainAll(system_soft,domain_decompose_weight);
1741  //std::cout<<"rank: "<<my_rank<<" weight: "<<domain_decompose_weight<<std::endl;
1742  }
1743 #ifdef PROFILE
1744  profile.domain.barrier();
1745  PS::Comm::barrier();
1746  profile.domain.end();
1747 #endif
1748  }
1749 
1751 
1754  bool adjustDtTreeReduce(PS::F64& _dt_reduce_factor, const PS::F64 _dt_tree_base, const PS::F64 _dt_tree_now) {
1755  bool dt_mod_flag = false;
1756  PS::F64 dt_tree_new = _dt_tree_base / _dt_reduce_factor;
1757  if(_dt_tree_now!= dt_tree_new) {
1758  // in increasing case, need to make sure the time is consistent with block step time/dt_tree
1759  if(_dt_tree_now<dt_tree_new) {
1760  PS::F64 dt_tree_max = calcDtLimit(stat.time, _dt_tree_base, _dt_tree_now);
1761  if (dt_tree_max == _dt_tree_now) {
1762  // no modificatoin, recover original reduce factor
1763  _dt_reduce_factor= _dt_tree_base/_dt_tree_now;
1764  }
1765  else {
1766  dt_mod_flag = true;
1767  // limit dt_reduce_factor to the maximum block step allown
1768  _dt_reduce_factor = std::min(_dt_tree_base/dt_tree_max, _dt_reduce_factor);
1769 
1770  }
1771  }
1772  else dt_mod_flag = true;
1773  }
1774 
1775  return dt_mod_flag;
1776  }
1777 
1779  void updateStatus(const bool _initial_flag) {
1780 #ifdef PROFILE
1781  profile.status.start();
1782 #endif
1783 // stat.shiftToOriginFrame(&system_soft[0], stat.n_real_loc);
1784  if (_initial_flag) {
1785  // calculate initial energy
1786  stat.energy.clear();
1787 #ifdef RECORD_CM_IN_HEADER
1788  stat.energy.calc(&system_soft[0], stat.n_real_loc, true, &(stat.pcm.pos), &(stat.pcm.vel));
1789 #else
1790  stat.energy.calc(&system_soft[0], stat.n_real_loc, true);
1791 #endif
1792 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1793  stat.energy.getSumMultiNodes(true);
1794 #endif
1795 #ifdef HARD_CHECK_ENERGY
1796  stat.energy.ekin_sd = stat.energy.ekin;
1797  stat.energy.epot_sd = stat.energy.epot;
1798 #endif
1799 #ifndef RECORD_CM_IN_HEADER
1800  stat.calcCenterOfMass(&system_soft[0], stat.n_real_loc);
1801 #endif
1802  //Ptcl::vel_cm = stat.pcm.vel;
1803 
1804  }
1805  else {
1806 #ifdef HARD_CHECK_ENERGY
1807  // reset sd energy reference to no slowdown case
1808  PS::F64 energy_old = stat.energy.ekin + stat.energy.epot;
1809  stat.energy.etot_sd_ref -= stat.energy.ekin_sd + stat.energy.epot_sd - energy_old;
1810 #endif
1811 #ifdef RECORD_CM_IN_HEADER
1812  stat.energy.calc(&system_soft[0], stat.n_real_loc,false, &(stat.pcm.pos), &(stat.pcm.vel));
1813 #else
1814  stat.energy.calc(&system_soft[0], stat.n_real_loc);
1815 #endif
1816 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1817  stat.energy.getSumMultiNodes();
1818 #endif
1819 
1820 #ifdef HARD_CHECK_ENERGY
1821  HardEnergy energy_local = system_hard_one_cluster.energy;
1822  energy_local += system_hard_isolated.energy;
1823 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1824  energy_local += system_hard_connected.energy;
1825 #endif
1826  // hard energy error
1827  stat.energy.error_hard_cum += PS::Comm::getSum(energy_local.de);
1828  stat.energy.error_hard_sd_cum += PS::Comm::getSum(energy_local.de_sd);
1829  // energy correction due to slowdown
1830  PS::F64 ekin_sd_correction = PS::Comm::getSum(energy_local.ekin_sd_correction);
1831  stat.energy.ekin_sd = stat.energy.ekin + ekin_sd_correction;
1832  PS::F64 epot_sd_correction = PS::Comm::getSum(energy_local.epot_sd_correction);
1833  stat.energy.epot_sd = stat.energy.epot + epot_sd_correction;
1834  PS::F64 etot_sd_correction = ekin_sd_correction + epot_sd_correction;
1835  PS::F64 de_change_cum = PS::Comm::getSum(energy_local.de_change_cum);
1836  PS::F64 de_change_binary_interrupt = PS::Comm::getSum(energy_local.de_change_binary_interrupt);
1837  PS::F64 de_change_modify_single = PS::Comm::getSum(energy_local.de_change_modify_single);
1838  stat.energy.etot_ref += de_change_cum;
1839  stat.energy.de_change_cum += de_change_cum;
1840  stat.energy.de_change_binary_interrupt += de_change_binary_interrupt;
1841  stat.energy.de_change_modify_single += de_change_modify_single;
1842  // for total energy reference, first add the cumulative change due to slowdown change in the integration (referring to no slowdown case), then add slowdown energy correction from current time
1843  PS::F64 de_sd_change_cum = PS::Comm::getSum(energy_local.de_sd_change_cum) + etot_sd_correction;
1844  PS::F64 de_sd_change_binary_interrupt = PS::Comm::getSum(energy_local.de_sd_change_binary_interrupt);
1845  PS::F64 de_sd_change_modify_single = PS::Comm::getSum(energy_local.de_sd_change_modify_single);
1846  stat.energy.etot_sd_ref += de_sd_change_cum;
1847  stat.energy.de_sd_change_cum += de_sd_change_cum;
1848  stat.energy.de_sd_change_binary_interrupt += de_sd_change_binary_interrupt;
1849  stat.energy.de_sd_change_modify_single += de_sd_change_modify_single;
1850 
1851  system_hard_one_cluster.energy.clear();
1852  system_hard_isolated.energy.clear();
1853 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1854  system_hard_connected.energy.clear();
1855 #endif
1856 #endif
1857 #ifndef RECORD_CM_IN_HEADER
1858  stat.calcCenterOfMass(&system_soft[0], stat.n_real_loc);
1859 #endif
1860  //Ptcl::vel_cm = stat.pcm.vel;
1861  }
1862  //stat.shiftToCenterOfMassFrame(&system_soft[0], stat.n_real_loc);
1863 
1864 #ifdef PROFILE
1865  profile.status.barrier();
1866  PS::Comm::barrier();
1867  profile.status.end();
1868 #endif
1869  }
1870 
1872  void output() {
1873 #ifdef PROFILE
1874  profile.output.start();
1875 #endif
1876  bool print_flag = input_parameters.print_flag;
1877  int write_style = input_parameters.write_style.value;
1878  std::cout<<std::setprecision(PRINT_PRECISION);
1879 
1880  // print status
1881  if(print_flag) {
1882  std::cout<<std::endl;
1883  stat.print(std::cout);
1884  }
1885 #ifdef GALPY
1886  if (print_flag) galpy_manager.printData(std::cout);
1887 #endif
1888  // write status, output to separate snapshots
1889  if(write_style==1) {
1890 
1891  // data output
1892  file_header.n_body = stat.n_real_glb;
1893  file_header.time = stat.time;
1894  file_header.nfile++;
1895 #ifdef RECORD_CM_IN_HEADER
1896  file_header.pos_offset = stat.pcm.pos;
1897  file_header.vel_offset = stat.pcm.vel;
1898 #endif
1899 
1900  std::string fname = input_parameters.fname_snp.value+"."+std::to_string(file_header.nfile);
1901 #ifdef PETAR_DEBUG
1902  assert(system_soft.getNumberOfParticleLocal()== stat.n_all_loc);
1903 #endif
1904  system_soft.setNumberOfParticleLocal(stat.n_real_loc);
1905  if (input_parameters.data_format.value==1||input_parameters.data_format.value==3)
1906  system_soft.writeParticleAscii(fname.c_str(), file_header);
1907  else if(input_parameters.data_format.value==0||input_parameters.data_format.value==2)
1908  system_soft.writeParticleBinary(fname.c_str(), file_header);
1909  system_soft.setNumberOfParticleLocal(stat.n_all_loc);
1910 
1911  if(my_rank==0) {
1912  // status output
1913  stat.printColumn(fstatus, WRITE_WIDTH);
1914  fstatus<<std::endl;
1915 
1916 #ifdef GALPY
1917  // for External potential
1918  galpy_manager.writePotentialPars(fname+".galpy", stat.time);
1919 #endif
1920  }
1921  }
1922  // write all information in to fstatus
1923  else if(write_style==2&&my_rank==0) {
1924  // write snapshot with one line
1925  stat.printColumn(fstatus, WRITE_WIDTH);
1926  for (int i=0; i<stat.n_real_loc; i++) {
1927 #ifdef RECORD_CM_IN_HEADER
1928  system_soft[i].printColumnWithOffset(stat.pcm, fstatus, WRITE_WIDTH);
1929 #else
1930  system_soft[i].printColumn(fstatus, WRITE_WIDTH);
1931 #endif
1932  }
1933  fstatus<<std::endl;
1934  }
1935  // write status only
1936  else if(write_style==3&&my_rank==0) {
1937  stat.printColumn(fstatus, WRITE_WIDTH);
1938  fstatus<<std::endl;
1939  }
1940 
1941  // save current error
1942  stat.energy.saveEnergyError();
1943 
1944 #ifdef PROFILE
1945  profile.output.barrier();
1946  PS::Comm::barrier();
1947  profile.output.end();
1948 #endif
1949  }
1950 
1951 #ifdef PROFILE
1952  void calcProfile() {
1954 
1955  // profile analysis
1956 
1957  PS::S32 n_hard_single = system_hard_one_cluster.getPtcl().size();
1958  PS::S32 n_hard_isolated = system_hard_isolated.getPtcl().size();
1959 
1960  n_count.hard_single += n_hard_single;
1961  n_count.hard_isolated += n_hard_isolated;
1962 
1963  n_count_sum.hard_single += PS::Comm::getSum(n_hard_single);
1964  n_count_sum.hard_isolated += PS::Comm::getSum(n_hard_isolated);
1965 
1966  PS::S64 ARC_substep_sum = system_hard_isolated.ARC_substep_sum;
1967  PS::S64 ARC_tsyn_step_sum = system_hard_isolated.ARC_tsyn_step_sum;
1968  PS::S64 ARC_n_groups = system_hard_isolated.ARC_n_groups;
1969  PS::S64 ARC_n_groups_iso = system_hard_isolated.ARC_n_groups_iso;
1970  PS::S64 H4_step_sum = system_hard_isolated.H4_step_sum;
1971 #ifdef HARD_COUNT_NO_NEIGHBOR
1972  PS::S64 n_neighbor_zero = system_hard_isolated.n_neighbor_zero;
1973 #endif
1974 
1975 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1976  PS::S32 n_hard_connected = system_hard_connected.getPtcl().size();
1977  n_count.hard_connected += n_hard_connected;
1978  n_count_sum.hard_connected += PS::Comm::getSum(n_hard_connected);
1979 
1980  ARC_substep_sum += system_hard_connected.ARC_substep_sum;
1981  ARC_tsyn_step_sum += system_hard_connected.ARC_tsyn_step_sum;
1982  ARC_n_groups += system_hard_connected.ARC_n_groups;
1983  ARC_n_groups_iso += system_hard_connected.ARC_n_groups_iso;
1984  H4_step_sum += system_hard_connected.H4_step_sum;
1985 #ifdef HARD_COUNT_NO_NEIGHBOR
1986  n_neighbor_zero+= system_hard_connected.n_neighbor_zero;
1987 #endif
1988 #endif
1989 
1990  n_count.ARC_substep_sum += ARC_substep_sum;
1991  n_count.ARC_tsyn_step_sum+= ARC_tsyn_step_sum;
1992  n_count.ARC_n_groups += ARC_n_groups;
1993  n_count.ARC_n_groups_iso += ARC_n_groups_iso;
1994  n_count.H4_step_sum += H4_step_sum;
1995 #ifdef HARD_COUNT_NO_NEIGHBOR
1996  n_count.n_neighbor_zero += n_neighbor_zero;
1997 #endif
1998 
1999  n_count_sum.ARC_substep_sum += PS::Comm::getSum(ARC_substep_sum);
2000  n_count_sum.ARC_tsyn_step_sum+= PS::Comm::getSum(ARC_tsyn_step_sum);
2001  n_count_sum.ARC_n_groups += PS::Comm::getSum(ARC_n_groups);
2002  n_count_sum.ARC_n_groups_iso += PS::Comm::getSum(ARC_n_groups_iso);
2003  n_count_sum.H4_step_sum += PS::Comm::getSum(H4_step_sum);
2004 #ifdef HARD_COUNT_NO_NEIGHBOR
2005  n_count_sum.n_neighbor_zero += PS::Comm::getSum(n_neighbor_zero);
2006 #endif
2007 
2008  system_hard_isolated.ARC_substep_sum = 0;
2009  system_hard_isolated.ARC_tsyn_step_sum=0;
2010  system_hard_isolated.ARC_n_groups = 0;
2011  system_hard_isolated.ARC_n_groups_iso = 0;
2012  system_hard_isolated.H4_step_sum = 0;
2013 #ifdef HARD_COUNT_NO_NEIGHBOR
2014  system_hard_isolated.n_neighbor_zero = 0;
2015 #endif
2016 
2017 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2018  system_hard_connected.ARC_substep_sum = 0;
2019  system_hard_connected.ARC_tsyn_step_sum=0;
2020  system_hard_connected.ARC_n_groups = 0;
2021  system_hard_connected.ARC_n_groups_iso = 0;
2022  system_hard_connected.H4_step_sum = 0;
2023 #ifdef HARD_COUNT_NO_NEIGHBOR
2024  system_hard_connected.n_neighbor_zero = 0;
2025 #endif
2026 #endif
2027 
2028  n_count.clearClusterCount();
2029  n_count.clusterCount(1, n_hard_single);
2030 
2031  const PS::S32 n_isolated_cluster = system_hard_isolated.getNumberOfClusters();
2032  n_count.cluster_isolated += n_isolated_cluster;
2033  n_count_sum.cluster_isolated += PS::Comm::getSum(n_isolated_cluster);
2034 
2035  const PS::S32* isolated_cluster_n_list = system_hard_isolated.getClusterNumberOfMemberList();
2036  for (PS::S32 i=0; i<n_isolated_cluster; i++) n_count.clusterCount(isolated_cluster_n_list[i]);
2037 
2038 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2039  const PS::S32 n_connected_cluster = system_hard_connected.getNumberOfClusters();
2040  n_count.cluster_connected += n_connected_cluster;
2041  n_count_sum.cluster_connected += PS::Comm::getSum(n_connected_cluster);
2042  const PS::S32* connected_cluster_n_list = system_hard_connected.getClusterNumberOfMemberList();
2043  for (PS::S32 i=0; i<n_connected_cluster; i++) n_count.clusterCount(connected_cluster_n_list[i]);
2044 #endif
2045 
2046 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2047  PS::S32 cluster_hist_size_loc = n_count.n_cluster.size();
2048  PS::S32 n_member_in_cluster_loc[cluster_hist_size_loc];
2049  PS::S32 n_cluster_loc[cluster_hist_size_loc];
2050  n_count.getherClusterCount(n_member_in_cluster_loc, n_cluster_loc, cluster_hist_size_loc);
2051  PS::S32 cluster_hist_size_recv[n_proc];
2052  PS::S32 cluster_hist_size_disp[n_proc+1];
2053  PS::Comm::allGather(&cluster_hist_size_loc, 1, cluster_hist_size_recv);
2054  cluster_hist_size_disp[0] = 0;
2055  for (PS::S32 i=0; i<n_proc; i++){
2056  cluster_hist_size_disp[i+1] = cluster_hist_size_recv[i] + cluster_hist_size_disp[i];
2057  }
2058  PS::S32 cluster_hist_size_total = cluster_hist_size_disp[n_proc];
2059  PS::S32 n_member_in_cluster_recv[cluster_hist_size_total];
2060  PS::S32 n_cluster_recv[cluster_hist_size_total];
2061  PS::Comm::allGatherV(n_member_in_cluster_loc, cluster_hist_size_loc,
2062  n_member_in_cluster_recv, cluster_hist_size_recv, cluster_hist_size_disp);
2063  PS::Comm::allGatherV(n_cluster_loc, cluster_hist_size_loc,
2064  n_cluster_recv, cluster_hist_size_recv, cluster_hist_size_disp);
2065  for (PS::S32 i=0; i<cluster_hist_size_total; i++) {
2066  n_count_sum.clusterCount(n_member_in_cluster_recv[i], n_cluster_recv[i]);
2067  }
2068 #else
2069  n_count_sum.addClusterCount(n_count);
2070 #endif
2071 
2072  dn_loop++;
2073 
2074  }
2075 
2077  void clearProfile() {
2078  profile.clear();
2079  tree_soft_profile.clear();
2080  tree_nb_profile.clear();
2081 #if defined(USE_GPU) && defined(GPU_PROFILE)
2082  gpu_profile.clear();
2083  gpu_counter.clear();
2084 #endif
2085  n_count.clear();
2086  n_count_sum.clear();
2087  dn_loop=0;
2088  }
2089 
2091  void printProfile() {
2092 
2093  const SysProfile& profile_min = profile.getMin();
2094 
2095  if(input_parameters.print_flag) {
2096  std::cout<<std::setprecision(5);
2097  std::cout<<"Tree step number: "<<dn_loop<<std::endl;
2098 
2099 
2100  std::cout<<"**** Wallclock time per step (local): [Min/Max]\n";
2101  //std::cout<<std::setw(PRINT_WIDTH)<<"Rank";
2102  profile.dumpName(std::cout);
2103  std::cout<<std::endl;
2104 
2105  //std::cout<<std::setw(PRINT_WIDTH)<<my_rank;
2106  profile_min.dump(std::cout,dn_loop);
2107  std::cout<<std::endl;
2108  profile.dump(std::cout,dn_loop);
2109  std::cout<<std::endl;
2110 
2111  std::cout<<"**** FDPS tree soft force time profile (local):\n";
2112  tree_soft_profile.dumpName(std::cout);
2113  std::cout<<std::endl;
2114  tree_soft_profile.dump(std::cout,dn_loop);
2115  std::cout<<std::endl;
2116 
2117  std::cout<<"**** Tree neighbor time profile (local):\n";
2118  tree_nb_profile.dumpName(std::cout);
2119  std::cout<<std::endl;
2120  tree_nb_profile.dump(std::cout,dn_loop);
2121  std::cout<<std::endl;
2122 
2123 #if defined(USE_GPU) && defined(GPU_PROFILE)
2124  std::cout<<"**** GPU time profile (local):\n";
2125  gpu_profile.dumpName(std::cout);
2126  gpu_counter.dumpName(std::cout);
2127  std::cout<<std::endl;
2128  gpu_profile.dump(std::cout,dn_loop);
2129  gpu_counter.dump(std::cout,dn_loop);
2130  std::cout<<std::endl;
2131 #endif
2132 
2133  std::cout<<"**** Number per step (global):\n";
2134  n_count_sum.dumpName(std::cout);
2135  std::cout<<std::endl;
2136  n_count_sum.dump(std::cout,dn_loop);
2137  std::cout<<std::endl;
2138 
2139  std::cout<<"**** Number of members in clusters (global):\n";
2140  n_count_sum.printHist(std::cout,dn_loop);
2141  }
2142 
2143  if(input_parameters.write_style.value>0) {
2144  fprofile<<std::setprecision(WRITE_PRECISION);
2145  fprofile<<std::setw(WRITE_WIDTH)<<my_rank;
2146  fprofile<<std::setw(WRITE_WIDTH)<<stat.time
2147  <<std::setw(WRITE_WIDTH)<<dn_loop
2148  <<std::setw(WRITE_WIDTH)<<stat.n_real_loc;
2149  profile.dump(fprofile, dn_loop, WRITE_WIDTH);
2150  profile.dumpBarrier(fprofile, dn_loop, WRITE_WIDTH);
2151  tree_soft_profile.dump(fprofile, dn_loop, WRITE_WIDTH);
2152  tree_nb_profile.dump(fprofile, dn_loop, WRITE_WIDTH);
2153 #if defined(USE_GPU) && defined(GPU_PROFILE)
2154  gpu_profile.dump(fprofile, dn_loop, WRITE_WIDTH);
2155  gpu_counter.dump(fprofile, dn_loop, WRITE_WIDTH);
2156 #endif
2157  n_count.dump(fprofile, dn_loop, WRITE_WIDTH);
2158  fprofile<<std::endl;
2159  }
2160  }
2161 #endif
2162 
2164  PS::S32 getParticleAdrFromID(const PS::S64 _id) {
2165  auto item = id_adr_map.find(_id);
2166  if (item==id_adr_map.end()) return -1;
2167  else return item->second;
2168  }
2169 
2171  void addParticleInIdAdrMap(FPSoft& _ptcl) {
2172  id_adr_map[_ptcl.id] = _ptcl.adr;
2173  }
2174 
2176  void reconstructIdAdrMap() {
2177  id_adr_map.clear();
2178  for (PS::S32 i=0; i<stat.n_real_loc; i++) id_adr_map[system_soft[i].id] = i;
2179  }
2180 
2181 #ifdef STELLAR_EVOLUTION
2182  void correctSoftPotMassChange() {
2184  // correct soft potential energy due to mass change
2185 #pragma omp parallel for
2186  for (int k=0; k<mass_modify_list.size(); k++) {
2187  PS::S32 i = mass_modify_list[k];
2188  auto& pi = system_soft[i];
2189  PS::F64 dpot = pi.dm*pi.pot_soft;
2190  stat.energy.etot_ref += dpot;
2191  stat.energy.de_change_cum += dpot;
2192  stat.energy.etot_sd_ref += dpot;
2193  stat.energy.de_sd_change_cum += dpot;
2194  pi.dm = 0.0;
2195  // ghost particle case, check in remove_particle instead
2196  //if(pi.mass==0.0&&pi.group_data.artificial.isUnused()) {
2197  // remove_list.push_back(i);
2198  //}
2199  }
2200  mass_modify_list.resizeNoInitialize(0);
2201  }
2202 #endif
2203 
2205  PS::S32 removeParticleFromIdAdrMap(const PS::S64 _id) {
2206  auto item = id_adr_map.find(_id);
2207  if (item==id_adr_map.end()) return -1;
2208  else {
2209  PS::S32 adr = item->second;
2210  id_adr_map.erase(item);
2211  return adr;
2212  }
2213  }
2214 
2216  void removeParticles() {
2218  assert(n_interrupt_glb==0);
2219 #ifdef PROFILE
2220  profile.other.start();
2221 #endif
2222 
2223  // first remove artificial particles
2224  system_soft.setNumberOfParticleLocal(stat.n_real_loc);
2225 
2226  // set flag for particles in remove_list
2227  PS::S32 n_remove = remove_list.size();
2228  for (PS::S32 i=0; i<n_remove; i++) {
2229  auto& pi = system_soft[remove_list[i]];
2230  pi.group_data.artificial.setParticleTypeToUnused(); // sign for removing
2231  // if mass is not zero, correct energy
2232  if (pi.mass>0) {
2233  PS::F64 dpot = pi.mass*pi.pot_tot;
2234  PS::F64 dkin = 0.5*pi.mass*(pi.vel*pi.vel);
2235  PS::F64 eloss = dpot + dkin;
2236  stat.energy.etot_ref -= eloss;
2237  stat.energy.de_change_cum -= eloss;
2238  stat.energy.etot_sd_ref -= eloss;
2239  stat.energy.de_sd_change_cum -= eloss;
2240  pi.mass = 0.0;
2241  }
2242  }
2243  remove_list.resizeNoInitialize(0);
2244 
2245  // check all particles escaper status and remove
2246  const PS::S32 num_thread = PS::Comm::getNumberOfThread();
2247  PS::ReallocatableArray<PS::S32> remove_list_thx[num_thread];
2248  for (PS::S32 i=0; i<num_thread; i++) remove_list_thx[i].resizeNoInitialize(0);
2249 
2250 #pragma omp parallel
2251  {
2252  const PS::S32 ith = PS::Comm::getThreadNum();
2253 #pragma omp for
2254  for (PS::S32 i=0; i<stat.n_real_loc; i++) {
2255  auto& pi = system_soft[i];
2256  if (escaper.isEscaper(pi,stat.pcm)) {
2257  remove_list_thx[ith].push_back(i);
2258  PS::F64 dpot = pi.mass*pi.pot_tot;
2259  PS::F64 dkin = 0.5*pi.mass*(pi.vel*pi.vel);
2260  PS::F64 eloss = dpot + dkin;
2261  stat.energy.etot_ref -= eloss;
2262  stat.energy.de_change_cum -= eloss;
2263  stat.energy.etot_sd_ref -= eloss;
2264  stat.energy.de_sd_change_cum -= eloss;
2265  }
2266  // Registered removed particles have already done energy correction
2267  else if (pi.mass==0.0&&pi.group_data.artificial.isUnused())
2268  remove_list_thx[ith].push_back(i);
2269  }
2270  }
2271 
2272  // gether remove indices
2273  int n_esc=0;
2274  for (PS::S32 i=0; i<num_thread; i++)
2275  for (PS::S32 k=0; k<remove_list_thx[i].size(); k++) {
2276  PS::S32 index=remove_list_thx[i][k];
2277  remove_list.push_back(index);
2278  remove_id_record.push_back(system_soft[index].id);
2279  if (system_soft[index].mass>0) {
2280  if (input_parameters.write_style.value>0) {
2281  fesc<<std::setw(WRITE_WIDTH)<<stat.time;
2282  system_soft[index].printColumn(fesc,WRITE_WIDTH);
2283  fesc<<std::endl;
2284  }
2285  n_esc++;
2286  }
2287  }
2288 
2289  // Remove particles
2290  n_remove = remove_list.size();
2291  system_soft.removeParticle(remove_list.getPointer(), remove_list.size());
2292 
2293  stat.n_escape_glb += PS::Comm::getSum(n_esc);
2294  stat.n_remove_glb += PS::Comm::getSum(n_remove);
2295 
2296  // reset particle number
2297  stat.n_real_loc = stat.n_real_loc-remove_list.size();
2298  system_soft.setNumberOfParticleLocal(stat.n_real_loc);
2299  stat.n_real_glb = system_soft.getNumberOfParticleGlobal();
2300  remove_list.resizeNoInitialize(0);
2301 
2302 #ifdef PETAR_DEBUG
2303 #pragma omp parallel for
2304  for(PS::S32 i=0; i<stat.n_real_loc; i++){
2305  assert(system_soft[i].mass!=0.0);
2306  assert(!std::isinf(system_soft[i].pos[0]));
2307  assert(!std::isnan(system_soft[i].pos[0]));
2308  assert(!std::isinf(system_soft[i].vel[0]));
2309  assert(!std::isnan(system_soft[i].vel[0]));
2310  }
2311 #endif
2312 
2313 #ifdef PROFILE
2314  profile.other.barrier();
2315  PS::Comm::barrier();
2316  profile.other.end();
2317 #endif
2318  }
2319 
2320 
2322  PS::S32 getNumberOfRemovedParticlesLocal() const {
2323  return remove_id_record.size();
2324  }
2325 
2327  PS::S32 getNumberOfRemovedParticlesGlobal() const {
2328  PS::S32 n_remove_loc = remove_id_record.size();
2329  return PS::Comm::getSum(n_remove_loc);
2330  }
2331 
2333  PS::S64* getRemovedIDListLocal() const {
2334  return remove_id_record.getPointer();
2335  }
2336 
2338  void clearRemovedIDList() {
2339  remove_id_record.resizeNoInitialize(0);
2340  }
2341 
2343  void exchangeParticle() {
2344  assert(n_interrupt_glb==0);
2345 
2346 #ifdef PROFILE
2347  profile.exchange.start();
2348 #endif
2349  system_soft.exchangeParticle(dinfo);
2350 
2351  const PS::S32 n_loc = system_soft.getNumberOfParticleLocal();
2352 
2353 #pragma omp parallel for
2354  for(PS::S32 i=0; i<n_loc; i++){
2355  system_soft[i].rank_org = my_rank;
2356  system_soft[i].adr = i;
2357  }
2358 
2359  // record real particle n_loc/glb
2360  stat.n_real_loc = n_loc;
2361 #ifdef PETAR_DEBUG
2362  assert(stat.n_real_glb==system_soft.getNumberOfParticleGlobal());
2363 #endif
2364 
2365 #ifdef PROFILE
2366  profile.exchange.barrier();
2367  PS::Comm::barrier();
2368  profile.exchange.end();
2369 #endif
2370  }
2371 
2373  static void initialFDPS(int argc, char *argv[]) {
2374  if (!initial_fdps_flag) {
2375  PS::Initialize(argc, argv);
2376  initial_fdps_flag = true;
2377  }
2378  }
2379 
2381  static void finalizeFDPS() {
2382  if (initial_fdps_flag) {
2383  PS::Finalize();
2384  initial_fdps_flag = false;
2385  }
2386  }
2387 
2388 #ifdef FDPS_COMM
2389  // create a MPI communicator
2395  void createCommunicator(const int n, const int rank[]) {
2396  assert(initial_fdps_flag);
2397  assert(!initial_parameters_flag);
2398  comm_info.create(n,rank);
2399  system_soft.setCommInfo(comm_info);
2400  dinfo.setCommInfo(comm_info);
2401  tree_nb.setCommInfo(comm_info);
2402  tree_soft.setCommInfo(comm_info);
2403  my_rank = comm_info.getRank();
2404  n_proc = comm_info.getNumberOfProc();
2405  }
2406 
2407  // set a MPI communicator
2412  void setCommunicator(const PS::CommInfo &_comm_info) {
2413  assert(initial_fdps_flag);
2414  assert(!initial_parameters_flag);
2415  comm_info = _comm_info;
2416  system_soft.setCommInfo(comm_info);
2417  dinfo.setCommInfo(comm_info);
2418  tree_nb.setCommInfo(comm_info);
2419  tree_soft.setCommInfo(comm_info);
2420  my_rank = comm_info.getRank();
2421  n_proc = comm_info.getNumberOfProc();
2422  }
2423 #endif
2424 
2426  void printLogo(std::ostream & fout) const {
2427  fout<<"\n ******************************************\n"
2428  <<" ██████╗ ███████╗████████╗ █████╗ ██████╗\n"
2429  <<" ██╔══██╗██╔════╝╚══██╔══╝██╔══██╗██╔══██╗\n"
2430  <<" ██████╔╝█████╗ ██║ ███████║██████╔╝\n"
2431  <<" ██╔═══╝ ██╔══╝ ██║ ██╔══██║██╔══██╗\n"
2432  <<" ██║ ███████╗ ██║ ██║ ██║██║ ██║\n"
2433  <<" ╚═╝ ╚══════╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝\n"
2434  <<" ******************************************\n"
2435  <<"Version: "<<GetVersion()
2436  <<std::endl;
2437  }
2438 
2440  void printReference(std::ostream & fout) const{
2441  fout<<"============== PeTar ================\n"
2442  <<" Online document: (preparing)\n"
2443  <<" GitHub page: https://github.com/lwang-astro/PeTar\n"
2444  <<" License: MIT\n"
2445  <<" Please cite the following papers when you use PeTar for any publications\n"
2446  <<" FDPS: see above FPDS Logo message\n";
2447  AR::printReference(fout);
2448  fout<<" PeTar: Wang L., Iwasawa M., Nitadori K., Makino J., 2020, MNRAS, 497, 536\n";
2449 #ifdef BSE_BASE
2450  BSEManager::printReference(fout);
2451 #endif
2452 #ifdef GALPY
2453  GalpyManager::printReference(fout);
2454 #endif
2455  fout<<" Copyright (C) 2017\n"
2456  <<" Long Wang, Masaki Iwasawa, Keigo Nitadori, Junichiro Makino and many others\n";
2457  fout<<"====================================="
2458  <<std::endl;
2459  }
2460 
2462  void printFeatures(std::ostream & fout) const {
2463 #ifdef USE_QUAD
2464  fout<<"Use quadrupole moment in tree force calculation\n";
2465 #endif
2466 
2467 #ifdef SOFT_PERT
2468 #ifdef TIDAL_TENSOR_3RD
2469  fout<<"Use 3rd order tidal tensor method\n";
2470 #else
2471  fout<<"Use 2nd order tidal tensor method\n";
2472 #endif
2473 #endif
2474 
2475 #ifdef ORBIT_SAMPLING
2476  fout<<"Use orbit-sampling method\n";
2477 #else
2478  fout<<"Use Pseudoparticle multipole method\n";
2479 #endif
2480 
2481 #ifdef STELLAR_EVOLUTION
2482  fout<<"Use stellar evolution method: ";
2483 #ifdef BSE_BASE
2484  fout<<BSEManager::getBSEName()<<std::endl;
2485 #else
2486  fout<<"Base\n";
2487 #endif
2488 #endif
2489 
2490 #ifdef GALPY
2491  fout<<"Use external potential: Galpy\n";
2492 #endif
2493 
2494 #ifdef KDKDK_2ND
2495  fout<<"Use 2nd-order KDKDK mode for tree step\n";
2496 #endif
2497 #ifdef KDKDK_4TH
2498  fout<<"Use 4th-order KDKDK mode for tree step\n";
2499 #endif
2500 
2501 #ifdef CLUSTER_VELOCITY
2502  fout<<"Use orbit-dependent neighbor criterion\n";
2503 #endif
2504 
2505 #ifdef HARD_CHECK_ENERGY
2506  fout<<"Check energy of short-range (hard) integration\n";
2507 #endif
2508 
2509 #ifdef HARD_COUNT_NO_NEIGHBOR
2510  fout<<"Count isolated particles in hard clusters\n";
2511 #endif
2512 
2513 #ifdef USE_SIMD
2514  fout<<"Use SIMD\n";
2515 #ifdef P3T_64BIT
2516  fout<<"Use 64 bit SIMD n";
2517 #endif
2518 #endif
2519 
2520 #ifdef USE_FUGAKU
2521  fout<<"Use Fugaku\n";
2522 #endif
2523 
2524 #ifdef USE_GPU
2525  fout<<"Use GPU\n";
2526 #endif
2527 
2528 #ifdef GPERF_PROFILE
2529  fout<<"Use gperftools\n";
2530 #endif
2531 
2532 #ifdef PROFILE
2533  fout<<"Calculate profile\n";
2534 #endif
2535 
2536 #ifdef GPU_PROFILE
2537  fout<<"Calculate GPU profile\n";
2538 #endif
2539 
2540  AR::printFeatures(fout);
2541  H4::printFeatures(fout);
2542  }
2543 
2545  void printDebugFeatures(std::ostream & fout) const {
2546 #ifdef PETAR_DEBUG
2547  fout<<"Debug mode: PETAR\n";
2548 #endif
2549 #ifdef STABLE_CHECK_DEBUG
2550  fout<<"Debug mode: STABLE_CHECK\n";
2551 #endif
2552 #ifdef CLUSTER_DEBUG
2553  fout<<"Debug mode: CLUSTER\n";
2554 #endif
2555 #ifdef ARTIFICIAL_PARTICLE_DEBUG
2556  fout<<"Debug mode: ARTIFICIAL_PARTICLE\n";
2557 #endif
2558 #ifdef HARD_DEBUG
2559  fout<<"Debug mode: HARD\n";
2560 #endif
2561 #ifdef NAN_CHECK_DEBUG
2562  fout<<"Debug mode: NAN_CHECK\n";
2563 #endif
2564  AR::printDebugFeatures(fout);
2565  H4::printDebugFeatures(fout);
2566  }
2567 
2569 
2574  int readParameters(int argc, char *argv[]) {
2575  // print reference
2576  if (my_rank==0) {
2577  printLogo(std::cerr);
2578  printReference(std::cerr);
2579  printFeatures(std::cerr);
2580  printDebugFeatures(std::cerr);
2581  std::cerr<<"====================================="<<std::endl;
2582  }
2583 
2584  //assert(initial_fdps_flag);
2585  assert(!read_parameters_flag);
2586  // reading parameters
2587  opterr = 0;
2588  read_parameters_flag = true;
2589  if (my_rank==0) input_parameters.print_flag=true;
2590  else input_parameters.print_flag=false;
2591  int read_flag = input_parameters.read(argc,argv);
2592 #ifdef BSE_BASE
2593  if (my_rank==0) bse_parameters.print_flag=true;
2594  else bse_parameters.print_flag=false;
2595  bse_parameters.read(argc,argv);
2596 #endif
2597 #ifdef GALPY
2598  if (my_rank==0) galpy_parameters.print_flag=true;
2599  else galpy_parameters.print_flag=false;
2600  galpy_parameters.read(argc,argv);
2601 #endif
2602 
2603  // help case, return directly
2604  if (read_flag==-1) {
2605  // avoid segmentation fault due to FDPS clear function bug
2606  tree_nb.initialize(input_parameters.n_glb.value, input_parameters.theta.value, input_parameters.n_leaf_limit.value, input_parameters.n_group_limit.value);
2607  tree_soft.initialize(input_parameters.n_glb.value, input_parameters.theta.value, input_parameters.n_leaf_limit.value, input_parameters.n_group_limit.value);
2608 
2609  return read_flag;
2610  }
2611 
2612  bool print_flag = input_parameters.print_flag;
2613  int write_style = input_parameters.write_style.value;
2614 
2615 #ifdef PROFILE
2616  if(print_flag) {
2617  std::cout<<"----- Parallelization information -----\n";
2618  std::cout<<"MPI processors: "<<n_proc<<std::endl;
2619  std::cout<<"OMP threads: "<<PS::Comm::getNumberOfThread()<<std::endl;
2620  }
2621 
2622  // open profile file
2623  if(write_style>0) {
2624  std::string fproname=input_parameters.fname_snp.value+".prof.rank."+std::to_string(my_rank);
2625  if(input_parameters.append_switcher.value==1) fprofile.open(fproname.c_str(),std::ofstream::out|std::ofstream::app);
2626  else {
2627  fprofile.open(fproname.c_str(),std::ofstream::out);
2628 
2629  fprofile<<std::setw(WRITE_WIDTH)<<"my_rank"
2630  <<std::setw(WRITE_WIDTH)<<"Time"
2631  <<std::setw(WRITE_WIDTH)<<"N_steps"
2632  <<std::setw(WRITE_WIDTH)<<"N_real_loc";
2633  profile.dumpName(fprofile, WRITE_WIDTH);
2634  profile.dumpBarrierName(fprofile, WRITE_WIDTH);
2635  tree_soft_profile.dumpName(fprofile, WRITE_WIDTH);
2636  tree_nb_profile.dumpName(fprofile, WRITE_WIDTH);
2637 #if defined(USE_GPU) && defined(GPU_PROFILE)
2638  gpu_profile.dumpName(fprofile, WRITE_WIDTH);
2639  gpu_counter.dumpName(fprofile, WRITE_WIDTH);
2640 #endif
2641  n_count.dumpName(fprofile, WRITE_WIDTH);
2642  fprofile<<std::endl;
2643  }
2644  fprofile<<std::setprecision(WRITE_PRECISION);
2645  }
2646 #endif
2647 
2648  // open output files
2649  // status information output
2650  std::string& fname_snp = input_parameters.fname_snp.value;
2651  if(write_style>0&&my_rank==0) {
2652  if(input_parameters.append_switcher.value==1)
2653  fstatus.open((fname_snp+".status").c_str(),std::ofstream::out|std::ofstream::app);
2654  else {
2655  fstatus.open((fname_snp+".status").c_str(),std::ofstream::out);
2656  // write titles of columns
2657  stat.printColumnTitle(fstatus,WRITE_WIDTH);
2658  if (write_style==2) {
2659  for (int i=0; i<stat.n_real_loc; i++) system_soft[0].printColumnTitle(fstatus, WRITE_WIDTH);
2660  }
2661  fstatus<<std::endl;
2662  }
2663  fstatus<<std::setprecision(WRITE_PRECISION);
2664  }
2665 
2666  if(write_style>0) {
2667  // open escaper file
2668  std::string my_rank_str = std::to_string(my_rank);
2669  std::string fname_esc = fname_snp + ".esc." + my_rank_str;
2670  if(input_parameters.append_switcher.value==1)
2671  fesc.open(fname_esc.c_str(), std::ofstream::out|std::ofstream::app);
2672  else {
2673  fesc.open(fname_esc.c_str(), std::ofstream::out);
2674  // write titles of columns, will cause issue when gether different MPI ranks
2675  // fesc<<std::setw(WRITE_WIDTH)<<"Time";
2676  // FPSoft::printColumnTitle(fesc,WRITE_WIDTH);
2677  // fesc<<std::endl;
2678  }
2679  fesc<<std::setprecision(WRITE_PRECISION);
2680 
2681 #ifdef BSE_BASE
2682  // open SSE/BSE file
2683  std::string fsse_name = fname_snp + fsse_par_suffix + "." + my_rank_str;
2684  std::string fbse_name = fname_snp + fbse_par_suffix + "." + my_rank_str;
2685  if(input_parameters.append_switcher.value==1) {
2686  hard_manager.ar_manager.interaction.fout_sse.open(fsse_name.c_str(), std::ofstream::out|std::ofstream::app);
2687  hard_manager.ar_manager.interaction.fout_bse.open(fbse_name.c_str(), std::ofstream::out|std::ofstream::app);
2688  }
2689  else {
2690  hard_manager.ar_manager.interaction.fout_sse.open(fsse_name.c_str(), std::ofstream::out);
2691  hard_manager.ar_manager.interaction.fout_bse.open(fbse_name.c_str(), std::ofstream::out);
2692  }
2693  hard_manager.ar_manager.interaction.fout_sse<<std::setprecision(WRITE_PRECISION);
2694  hard_manager.ar_manager.interaction.fout_bse<<std::setprecision(WRITE_PRECISION);
2695 #endif
2696 
2697 #ifdef ADJUST_GROUP_PRINT
2698  // open file for new/end group information
2699  if (input_parameters.adjust_group_write_option.value==1) {
2700  std::string fgroup_name = fname_snp + ".group." + my_rank_str;
2701  if(input_parameters.append_switcher.value==1)
2702  hard_manager.h4_manager.fgroup.open(fgroup_name.c_str(), std::ofstream::out|std::ofstream::app);
2703  else
2704  hard_manager.h4_manager.fgroup.open(fgroup_name.c_str(), std::ofstream::out);
2705  hard_manager.h4_manager.fgroup<<std::setprecision(WRITE_PRECISION);
2706  }
2707 #endif
2708  }
2709 
2710 #ifdef HARD_DUMP
2711  // initial hard_dump
2712  const PS::S32 num_thread = PS::Comm::getNumberOfThread();
2713  hard_dump.initial(num_thread, my_rank);
2714 #endif
2715 
2716  // particle system
2717  system_soft.initialize();
2718 
2719  // domain decomposition
2720  system_soft.setAverageTargetNumberOfSampleParticlePerProcess(input_parameters.n_smp_ave.value);
2721  const PS::F32 coef_ema = 0.2;
2722  domain_decompose_weight=1.0;
2723  dinfo.initialize(coef_ema);
2724 
2725  if(pos_domain==NULL) {
2726  pos_domain = new PS::F64ort[n_proc];
2727  for(PS::S32 i=0; i<n_proc; i++) pos_domain[i] = dinfo.getPosDomain(i);
2728  }
2729 
2730  // tree for neighbor search
2731  tree_nb.initialize(input_parameters.n_glb.value, input_parameters.theta.value, input_parameters.n_leaf_limit.value, input_parameters.n_group_limit.value);
2732 
2733  // tree for force
2734  PS::S64 n_tree_init = input_parameters.n_glb.value + input_parameters.n_bin.value;
2735  tree_soft.initialize(n_tree_init, input_parameters.theta.value, input_parameters.n_leaf_limit.value, input_parameters.n_group_limit.value);
2736 
2737  // initial search cluster
2738  search_cluster.initialize();
2739 
2740  return read_flag;
2741  }
2742 
2744  void readDataFromFile() {
2745  assert(read_parameters_flag);
2746  if(input_parameters.fname_inp.value=="__NONE__") {
2747  std::cerr<<"Error: the input data filename is not provided, make sure your options have the correct format."<<std::endl;
2748  abort();
2749  }
2750 
2751  PS::S32 data_format = input_parameters.data_format.value;
2752  auto* data_filename = input_parameters.fname_inp.value.c_str();
2753 
2754  if(data_format==1||data_format==2||data_format==4)
2755  system_soft.readParticleAscii(data_filename, file_header);
2756  else
2757  system_soft.readParticleBinary(data_filename, file_header);
2758  PS::Comm::broadcast(&file_header, 1, 0);
2759  PS::S64 n_glb = system_soft.getNumberOfParticleGlobal();
2760  PS::S64 n_loc = system_soft.getNumberOfParticleLocal();
2761 
2762  if(input_parameters.print_flag) {
2763  std::cout<<std::setprecision(WRITE_PRECISION);
2764  std::cout<<"----- Reading file: "<<data_filename<<" -----"<<std::endl
2765  <<"Number of particles = "<<n_glb<<std::endl
2766  <<"Time = "<<file_header.time<<std::endl;
2767  }
2768 
2769  stat.time = file_header.time;
2770  stat.n_real_glb = stat.n_all_glb = n_glb;
2771  stat.n_real_loc = stat.n_all_loc = n_loc;
2772 #ifdef RECORD_CM_IN_HEADER
2773  PS::F64 mass = 0.0;
2774  for (int i=0; i<n_loc; i++) mass += system_soft[i].mass;
2775 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2776  stat.pcm.mass = PS::Comm::getSum(mass);
2777 #else
2778  stat.pcm.mass = mass;
2779 #endif
2780  stat.pcm.pos = file_header.pos_offset;
2781  stat.pcm.vel = file_header.vel_offset;
2782  stat.pcm.is_center_shift_flag = true;
2783 #endif
2784 
2785  input_parameters.n_glb.value = n_glb;
2786 
2787  read_data_flag = true;
2788  }
2789 
2791 
2795  template <class Tptcl>
2796  void readDataFromArray(PS::S64 _n_partcle, Tptcl* _particle) {
2797  assert(read_parameters_flag);
2798 
2799  PS::S64 n_glb = _n_partcle;
2800 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2801  PS::S64 n_loc = n_glb / n_proc;
2802  if( n_glb % n_proc > my_rank) n_loc++;
2803  PS::S64 i_h = n_glb/n_proc*my_rank;
2804  if( n_glb % n_proc > my_rank) i_h += my_rank;
2805  else i_h += n_glb % n_proc;
2806 #else
2807  PS::S64 n_loc = n_glb;
2808  PS::S64 i_h = 0;
2809 #endif
2810 
2811  system_soft.setNumberOfParticleLocal(n_loc);
2812  for(PS::S32 i=0; i<n_loc; i++){
2813  int k = i_h+i;
2814  system_soft[i].mass = _particle[k].getMass();
2815  system_soft[i].pos.x = _particle[k].getPos()[0];
2816  system_soft[i].pos.y = _particle[k].getPos()[1];
2817  system_soft[i].pos.z = _particle[k].getPos()[2];
2818  system_soft[i].vel.x = _particle[k].getVel()[0];
2819  system_soft[i].vel.y = _particle[k].getVel()[1];
2820  system_soft[i].vel.z = _particle[k].getVel()[2];
2821  system_soft[i].id = i_h+i+1;
2822  system_soft[i].group_data.artificial.setParticleTypeToSingle();
2823  }
2824 
2825  file_header.nfile = 0;
2826  file_header.n_body = n_glb;
2827  file_header.time = 0.0;
2828 
2829  stat.time = 0.0;
2830  stat.n_real_glb = stat.n_all_glb = n_glb;
2831  stat.n_real_loc = stat.n_all_loc = n_loc;
2832 #ifdef RECORD_CM_IN_HEADER
2833  stat.calcAndShiftCenterOfMass(&system_soft[0], n_loc, 3, true);
2834  file_header.pos_offset = stat.pcm.pos;
2835  file_header.vel_offset = stat.pcm.vel;
2836 #endif
2837  input_parameters.n_glb.value = n_glb;
2838 
2839  read_data_flag = true;
2840  }
2841 
2843 
2854  void readDataFromArray(PS::S64 _n_partcle, PS::F64* _mass, PS::F64* _x, PS::F64* _y, PS::F64* _z, PS::F64* _vx, PS::F64* _vy, PS::F64* _vz, PS::S64* _id=NULL) {
2855  assert(read_parameters_flag);
2856 
2857  PS::S64 n_glb = _n_partcle;
2858 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
2859  PS::S64 n_loc = n_glb / n_proc;
2860  if( n_glb % n_proc > my_rank) n_loc++;
2861  PS::S64 i_h = n_glb/n_proc*my_rank;
2862  if( n_glb % n_proc > my_rank) i_h += my_rank;
2863  else i_h += n_glb % n_proc;
2864 #else
2865  PS::S64 n_loc = n_glb;
2866  PS::S64 i_h = 0;
2867 #endif
2868 
2869  system_soft.setNumberOfParticleLocal(n_loc);
2870  if (_id!=NULL) {
2871  for(PS::S32 i=0; i<n_loc; i++){
2872  system_soft[i].mass = _mass[i_h+i];
2873  system_soft[i].pos.x = _x[i_h+i];
2874  system_soft[i].pos.y = _y[i_h+i];
2875  system_soft[i].pos.z = _z[i_h+i];
2876  system_soft[i].vel.x = _vx[i_h+i];
2877  system_soft[i].vel.y = _vy[i_h+i];
2878  system_soft[i].vel.z = _vz[i_h+i];
2879  system_soft[i].id = _id[i_h+i];
2880  system_soft[i].group_data.artificial.setParticleTypeToSingle();
2881  assert(_id[i_h+i]>0);
2882  }
2883  }
2884  else {
2885  for(PS::S32 i=0; i<n_loc; i++){
2886  system_soft[i].mass = _mass[i_h+i];
2887  system_soft[i].pos.x = _x[i_h+i];
2888  system_soft[i].pos.y = _y[i_h+i];
2889  system_soft[i].pos.z = _z[i_h+i];
2890  system_soft[i].vel.x = _vx[i_h+i];
2891  system_soft[i].vel.y = _vy[i_h+i];
2892  system_soft[i].vel.z = _vz[i_h+i];
2893  system_soft[i].id = i_h+i+1;
2894  system_soft[i].group_data.artificial.setParticleTypeToSingle();
2895  }
2896  }
2897  file_header.nfile = 0;
2898  file_header.n_body = n_glb;
2899  file_header.time = 0.0;
2900 
2901  stat.time = 0.0;
2902  stat.n_real_glb = stat.n_all_glb = n_glb;
2903  stat.n_real_loc = stat.n_all_loc = n_loc;
2904 #ifdef RECORD_CM_IN_HEADER
2905  stat.calcAndShiftCenterOfMass(&system_soft[0], n_loc, 3, true);
2906 
2907  file_header.pos_offset = stat.pcm.pos;
2908  file_header.vel_offset = stat.pcm.vel;
2909 #endif
2910 
2911  input_parameters.n_glb.value = n_glb;
2912 
2913  read_data_flag = true;
2914  }
2915 
2917  void generatePlummer() {
2918  // ensure parameters are used
2919  assert(read_parameters_flag);
2920 
2921  PS::S64 n_glb = input_parameters.n_glb.value;
2922  assert(n_glb>0);
2923 
2924  PS::S64 n_loc = n_glb / n_proc;
2925  if( n_glb % n_proc > my_rank) n_loc++;
2926  system_soft.setNumberOfParticleLocal(n_loc);
2927 
2928  PS::F64 * mass;
2929  PS::F64vec * pos;
2930  PS::F64vec * vel;
2931 
2932  const PS::F64 m_tot = 1.0;
2933  const PS::F64 eng = -0.25;
2934 
2935  ParticleDistributionGenerator::makePlummerModel(m_tot, n_glb, n_loc, mass, pos, vel, eng);
2936 
2937  PS::S64 i_h = n_glb/n_proc*my_rank;
2938  if( n_glb % n_proc > my_rank) i_h += my_rank;
2939  else i_h += n_glb % n_proc;
2940  for(PS::S32 i=0; i<n_loc; i++){
2941  system_soft[i].mass = mass[i];
2942  system_soft[i].pos = pos[i];
2943  system_soft[i].vel = vel[i];
2944  system_soft[i].id = i_h + i + 1;
2945 #ifdef STELLAR_EVOLUTION
2946  system_soft[i].radius = 0.0;
2947  system_soft[i].dm = 0.0;
2948  system_soft[i].time_record = 0.0;
2949  system_soft[i].time_interrupt = 0.0;
2950  system_soft[i].binary_state = 0;
2951 #ifdef BSE_BASE
2952  system_soft[i].star.initial(mass[i]*bse_parameters.mscale.value);
2953 #endif
2954 #endif
2955  system_soft[i].group_data.artificial.setParticleTypeToSingle();
2956  }
2957 
2958  file_header.nfile = 0;
2959  file_header.n_body = n_glb;
2960  file_header.time = 0.0;
2961 
2962  stat.time = 0.0;
2963  stat.n_real_glb = stat.n_all_glb = n_glb;
2964  stat.n_real_loc = stat.n_all_loc = n_loc;
2965 
2966 #ifdef RECORD_CM_IN_HEADER
2967  stat.calcAndShiftCenterOfMass(&system_soft[0], n_loc, 1, true);
2968  file_header.pos_offset = stat.pcm.pos;
2969  file_header.vel_offset = stat.pcm.vel;
2970 #endif
2971 
2972  read_data_flag = true;
2973 
2974  delete [] mass;
2975  delete [] pos;
2976  delete [] vel;
2977  }
2978 
2980 
2982  void generateKeplerDisk(const PS::F64 ax_in, // [AU]
2983  const PS::F64 ax_out, // [AU]
2984  const PS::F64 ecc_rms, // normalized
2985  const PS::F64 inc_rms, // normalized
2986  const PS::F64 dens = 10.0, // [g/cm^2]
2987  const PS::F64 mass_sun = 1.0, //[m_sun]
2988  const double a_ice = 0.0,
2989  const double f_ice = 1.0,
2990  const double power = -1.5,
2991  const PS::S32 seed = 0) {
2992 
2993  // ensure parameters are used
2994  assert(read_parameters_flag);
2995 
2996  PS::S64 n_glb = input_parameters.n_glb.value;
2997  assert(n_glb>0);
2998 
2999  PS::S64 n_loc = n_glb / n_proc;
3000  if( n_glb % n_proc > my_rank) n_loc++;
3001  system_soft.setNumberOfParticleLocal(n_loc);
3002 
3003  PS::F64 * mass;
3004  PS::F64vec * pos;
3005  PS::F64vec * vel;
3006 
3007  PS::F64 mass_planet_glb;
3008  ParticleDistributionGenerator::makeKeplerDisk(mass_planet_glb, mass, pos, vel, n_glb, n_loc,
3009  ax_in, ax_out, ecc_rms, inc_rms, dens, mass_sun, a_ice, f_ice, power, seed);
3010 
3011  PS::S64 i_h = n_glb/n_proc*my_rank;
3012  if( n_glb % n_proc > my_rank) i_h += my_rank;
3013  else i_h += n_glb % n_proc;
3014  for(PS::S32 i=0; i<n_loc; i++){
3015  system_soft[i].mass = mass[i];
3016  system_soft[i].pos = pos[i];
3017  system_soft[i].vel = vel[i];
3018  system_soft[i].id = i_h + i + 1;
3019  system_soft[i].group_data.artificial.setParticleTypeToSingle();
3020  }
3021 
3022  file_header.nfile = 0;
3023  stat.time = file_header.time = 0.0;
3024  file_header.n_body = n_glb;
3025 
3026  stat.n_real_glb = stat.n_all_glb = n_glb;
3027  stat.n_real_loc = stat.n_all_loc = n_loc;
3028 
3029 #ifdef RECORD_CM_IN_HEADER
3030  stat.calcAndShiftCenterOfMass(&system_soft[0], n_loc, 1, true);
3031  file_header.pos_offset = stat.pcm.pos;
3032  file_header.vel_offset = stat.pcm.vel;
3033 #endif
3034 
3035  read_data_flag = true;
3036 
3037  delete [] mass;
3038  delete [] pos;
3039  delete [] vel;
3040  }
3041 
3043  void initialParameters() {
3044  // ensure data is read
3045  assert(read_data_flag);
3046  assert(stat.n_real_glb>0);
3047 
3048  bool print_flag = input_parameters.print_flag;
3049  int write_style = input_parameters.write_style.value;
3050  std::cout<<std::setprecision(WRITE_PRECISION);
3051 
3052  // units
3053  if (input_parameters.unit_set.value==1) {
3054  input_parameters.gravitational_constant.value = 0.00449830997959438; // pc^3/(Msun*Myr^2)
3055 #ifdef BSE_BASE
3056  bse_parameters.tscale.value = 1.0; // Myr
3057  bse_parameters.rscale.value = 44353565.919218; // pc -> rsun
3058  bse_parameters.mscale.value = 1.0; // Msun
3059  bse_parameters.vscale.value = 0.977813107686401; // pc/Myr -> km/s
3060 #endif
3061  if(print_flag) {
3062  std::cout<<"----- Unit set 1: Msun, pc, Myr -----\n"
3063  <<"gravitational_constant = "<<input_parameters.gravitational_constant.value<<" pc^3/(Msun*Myr^2)\n";
3064 #ifdef BSE_BASE
3065  std::cout<<"----- Unit conversion for BSE ----- \n"
3066  <<" tscale = "<<bse_parameters.tscale.value<<" Myr / Myr\n"
3067  <<" mscale = "<<bse_parameters.mscale.value<<" Msun / Msun\n"
3068  <<" rscale = "<<bse_parameters.rscale.value<<" Rsun / pc\n"
3069  <<" vscale = "<<bse_parameters.vscale.value<<" [km/s] / [pc/Myr]\n";
3070 #endif
3071 
3072  }
3073 //#ifdef GALPY
3074 // galpy_parameters.setStdUnit(print_flag);
3075 //#endif
3076 
3077  }
3078 
3079  // calculate system parameters
3080  PS::F64 r_in, mass_average, vel_disp;// mass_max, vel_max;
3081  PS::F64& r_out = input_parameters.r_out.value;
3082  PS::F64& r_bin = input_parameters.r_bin.value;
3083  PS::F64& r_search_min = input_parameters.r_search_min.value;
3084 // PS::F64& r_search_max = input_parameters.r_search_max.value;
3085  PS::F64& dt_soft = input_parameters.dt_soft.value;
3086  PS::F64& dt_snap = input_parameters.dt_snap.value;
3087  PS::F64& search_vel_factor = input_parameters.search_vel_factor.value;
3088  PS::F64& ratio_r_cut = input_parameters.ratio_r_cut.value;
3089  PS::S64& n_bin = input_parameters.n_bin.value;
3090  //PS::F64& theta = input_parameters.theta.value;
3091  PS::F64& G = input_parameters.gravitational_constant.value;
3092  PS::F64& nstep_dt_soft_kepler = input_parameters.nstep_dt_soft_kepler.value;
3093 
3094  // local particle number
3095  const PS::S64 n_loc = system_soft.getNumberOfParticleLocal();
3096 
3097  // local c.m velocity
3098  PS::F64vec vel_cm_loc = 0.0;
3099  // local c.m. pos
3100  PS::F64vec pos_cm_loc = 0.0;
3101  // local c.m. mass
3102  PS::F64 mass_cm_loc = 0.0;
3103  // local maximum mass
3104  //PS::F64 mass_max_loc = 0.0;
3105  // box size
3106  PS::F64 rmax=0.0;
3107 
3108  for(PS::S64 i=0; i<n_loc; i++){
3109  PS::F64 mi = system_soft[i].mass;
3110  PS::F64vec vi = system_soft[i].vel;
3111  PS::F64vec ri = system_soft[i].pos;
3112 
3113 #ifdef PETAR_DEBUG
3114  assert(mi>0);
3115 #endif
3116  mass_cm_loc += mi;
3117  vel_cm_loc += mi * vi;
3118  pos_cm_loc += mi * ri;
3119  PS::F64 r2 = ri*ri;
3120  rmax = std::max(r2,rmax);
3121  //mass_max_loc = std::max(mi, mass_max_loc);
3122  }
3123  rmax = std::sqrt(rmax);
3124  PS::F64 rmax_glb = PS::Comm::getMaxValue(rmax);
3125 
3126  // global c.m. parameters
3127  PS::F64 mass_cm_glb = PS::Comm::getSum(mass_cm_loc);
3128  //mass_max = PS::Comm::getMaxValue(mass_max_loc);
3129  PS::F64vec pos_cm_glb = PS::Comm::getSum(pos_cm_loc);
3130  PS::F64vec vel_cm_glb = PS::Comm::getSum(vel_cm_loc);
3131  pos_cm_glb /= mass_cm_glb;
3132  vel_cm_glb /= mass_cm_glb;
3133 
3134  PS::F64 rmin_glb = std::sqrt(pos_cm_glb*pos_cm_glb);
3135 
3136  // local velocity square
3137  PS::F64 vel_sq_loc = 0.0;
3138  PS::S64 n_vel_loc_count = 0;
3139 
3140  // single particle starting index
3141  PS::S64 single_start_index = 0;
3142  const PS::S64 bin_last_id = 2*n_bin;
3143  if (system_soft[0].id<bin_last_id) {
3144  single_start_index = std::min(bin_last_id - system_soft[0].id + 1,n_loc);
3145  if(single_start_index%2!=0) single_start_index--;
3146  }
3147  // binary particle starting index
3148  const PS::S64 binary_start_index = (system_soft[0].id%2==0)?1:0;
3149 
3150  // calculate velocity dispersion
3151  for (PS::S64 i=binary_start_index; i<single_start_index; i+=2) {
3152  PS::F64 m1 = system_soft[i].mass;
3153  PS::F64 m2 = system_soft[i+1].mass;
3154  PS::F64vec dv = (m1*system_soft[i].vel + m2*system_soft[i+1].vel)/(m1+m2) - vel_cm_glb;
3155  vel_sq_loc += dv * dv;
3156  n_vel_loc_count++;
3157  }
3158 
3159  if (single_start_index <n_loc)
3160  for (PS::S64 i=single_start_index; i<n_loc; i++){
3161  PS::F64vec dv = system_soft[i].vel - vel_cm_glb;
3162  vel_sq_loc += dv * dv;
3163  n_vel_loc_count++;
3164  }
3165 
3166  const PS::S64 n_vel_glb_count= PS::Comm::getSum(n_vel_loc_count);
3167  const PS::S64 n_glb = PS::Comm::getSum(n_loc);
3168  const PS::F64 vel_sq_glb = PS::Comm::getSum(vel_sq_loc);
3169  vel_disp = sqrt(vel_sq_glb / 3.0 / (PS::F64)n_vel_glb_count);
3170 
3171  PS::F64 mass_average_glb = mass_cm_glb/(PS::F64)n_glb;
3172  mass_average = mass_average_glb;
3173 
3174  // flag to check whether r_ous is already defined
3175  bool r_out_flag = (r_out>0);
3176 
3177  // if r_out is already defined, calculate r_in based on ratio_r_cut
3178  if (r_out_flag) r_in = r_out * ratio_r_cut;
3179  // calculate r_out based on virial radius scaled with (N)^(1/3), calculate r_in by ratio_r_cut
3180  else {
3181  if (n_glb>1) {
3182  r_out = std::min(0.1*G*mass_cm_glb/(std::pow(n_glb,1.0/3.0)) / (3*vel_disp*vel_disp), 3.0*(rmax_glb-rmin_glb));
3183  r_in = r_out * ratio_r_cut;
3184  }
3185  else {
3186  // give two small values, no meaning at all
3187  r_out = 1e-16;
3188  r_in = r_out*ratio_r_cut;
3189  if (print_flag) std::cout<<"In one particle case, changeover radius is set to a small value\n";
3190  }
3191  }
3192 
3193  // if tree time step is not defined, calculate tree time step by r_out and velocity dispersion
3194  if (dt_soft==0.0) {
3195  if (n_glb==1) {
3196  if (print_flag) std::cout<<"In one particle case, tree time step is finishing - starting time\n";
3197  dt_soft = input_parameters.time_end.value - stat.time;
3198  }
3199  else {
3200  // 1/nstep of a binary period with semi-major axis = r_int.
3201  if (nstep_dt_soft_kepler>0)
3202  dt_soft = regularTimeStep(COMM::Binary::semiToPeriod(r_in, mass_average, G)/nstep_dt_soft_kepler);
3203  else
3204  dt_soft = regularTimeStep(0.1*r_out / vel_disp);
3205  }
3206  }
3207  else {
3208  dt_soft = regularTimeStep(dt_soft);
3209  // if r_out is not defined, adjust r_out to minimum based on tree step
3210  if (!r_out_flag) {
3211  if (n_glb>1) {
3212  if (nstep_dt_soft_kepler>0) {
3213  r_in = COMM::Binary::periodToSemi(dt_soft*nstep_dt_soft_kepler, mass_average, G);
3214  r_out = r_in / ratio_r_cut;
3215  }
3216  else {
3217  r_out = 10.0*dt_soft*vel_disp;
3218  r_in = r_out * ratio_r_cut;
3219  }
3220  }
3221  else {
3222  r_out = 1e-16;
3223  r_in = r_out*ratio_r_cut;
3224  if (print_flag) std::cout<<"In one particle case, changeover radius is set to a small value\n";
3225  }
3226  }
3227  }
3228 
3229  // if r_bin is not defined, set to theta * r_in
3230  if (r_bin==0.0) r_bin = 0.8*r_in;
3231 
3232  // if r_search_min is not defined, calculate by search_vel_factor*velocity_dispersion*tree_time_step + r_out
3233  if (r_search_min==0.0) r_search_min = search_vel_factor*vel_disp*dt_soft + r_out;
3234  // if r_search_max is not defined, calcualte by 5*r_out
3235 // if (r_search_max==0.0) r_search_max = 5*r_out;
3236  // calculate v_max based on r_search_max, tree time step and search_vel_factor
3237  //vel_max = (r_search_max - r_out) / dt_soft / search_vel_factor;
3238 
3239  // regularize output time
3240  dt_snap = regularTimeStep(dt_snap);
3241 
3242  EPISoft::eps = input_parameters.eps.value;
3243  EPISoft::r_out = r_out;
3244  ForceSoft::grav_const = input_parameters.gravitational_constant.value;
3245  Ptcl::search_factor = search_vel_factor;
3246  Ptcl::r_search_min = r_search_min;
3247  Ptcl::mean_mass_inv = 1.0/mass_average;
3248  Ptcl::r_group_crit_ratio = r_bin/r_in;
3249  //Ptcl::vel_cm = vel_cm_glb;
3250  escaper.r_escape_sq = input_parameters.r_escape.value*input_parameters.r_escape.value;
3251  escaper.check_energy_flag = (input_parameters.r_escape.value>=0);
3252 
3253  if(print_flag) {
3254  // set print format
3255  std::cout<<"----- Parameter list: -----\n";
3256  std::cout<<" mass_average = "<<mass_average <<std::endl
3257  <<" r_in = "<<r_in <<std::endl
3258  <<" r_out = "<<r_out <<std::endl
3259  <<" r_bin = "<<r_bin <<std::endl
3260  <<" r_search_min = "<<r_search_min <<std::endl
3261  <<" vel_disp = "<<vel_disp <<std::endl
3262  <<" dt_soft = "<<dt_soft <<std::endl;
3263  }
3264 
3265  // check restart
3266  bool restart_flag = file_header.nfile; // nfile = 0 is assumed as initial data file
3267 
3268  // set id_offset
3269 #ifdef PETAR_DEBUG
3270  assert(stat.n_real_glb == input_parameters.n_glb.value);
3271  assert(stat.n_real_glb == system_soft.getNumberOfParticleGlobal());
3272 #endif
3273  PS::S64& id_offset = input_parameters.id_offset.value;
3274  id_offset = id_offset==-1 ? stat.n_real_glb+1 : id_offset;
3275 
3276  // initial particles paramters
3277  if (!restart_flag) {
3278 #pragma omp parallel for
3279  for (PS::S32 i=0; i<stat.n_real_loc; i++) {
3280  // ID safety check
3281  PS::S64 id = system_soft[i].id;
3282  assert(id_offset>id);
3283  assert(id>0);
3284 
3285  // set changeover
3286  PS::F64 m_fac = system_soft[i].mass*Ptcl::mean_mass_inv;
3287  system_soft[i].changeover.setR(m_fac, r_in, r_out);
3288 
3289  // calculate r_search for particles, for binary, r_search depend on vel_disp
3290  if(id<=2*n_bin) system_soft[i].r_search = std::max(r_search_min,vel_disp*dt_soft*search_vel_factor + system_soft[i].changeover.getRout());
3291  else system_soft[i].calcRSearch(dt_soft);
3292 
3293  auto& pi_cm = system_soft[i].group_data.cm;
3294  pi_cm.mass = pi_cm.vel.x = pi_cm.vel.y = pi_cm.vel.z = 0.0;
3295  }
3296  }
3297  else {
3298  // clear up group_data.cm to avoid issue in search neighbor; update changeover and r_search
3299 #pragma omp parallel for
3300  for (PS::S32 i=0; i<stat.n_real_loc; i++) {
3301  auto& pi_cm = system_soft[i].group_data.cm;
3302 
3303  // update changeover
3304  if (input_parameters.update_changeover_flag) {
3305  PS::F64 m_fac = system_soft[i].mass*Ptcl::mean_mass_inv;
3306  system_soft[i].changeover.setR(m_fac, r_in, r_out);
3307  }
3308 
3309  // update research
3310  if (input_parameters.update_rsearch_flag) {
3311  if (pi_cm.mass!=0.0) {
3312 //#ifdef GROUP_DATA_WRITE_ARTIFICIAL
3313  system_soft[i].r_search = std::max(r_search_min, vel_disp*dt_soft*search_vel_factor + system_soft[i].changeover.getRout());
3314 // not correct, when data is dumped, it is not cm information but artificial particle data
3315 //#else
3316 // PS::F64 vcm = std::sqrt(pi_cm.vel*pi_cm.vel);
3317 // system_soft[i].r_search = std::max(r_search_min,vcm*dt_soft*search_vel_factor + system_soft[i].changeover.getRout());
3318 //#endif
3319  }
3320  else system_soft[i].calcRSearch(dt_soft);
3321  }
3322 
3323  pi_cm.mass = pi_cm.vel.x = pi_cm.vel.y = pi_cm.vel.z = 0.0;
3324  }
3325 #ifdef RECORD_CM_IN_HEADER
3326  stat.calcAndShiftCenterOfMass(&system_soft[0], stat.n_real_loc);
3327 #else
3328  stat.calcCenterOfMass(&system_soft[0], stat.n_real_loc);
3329 #endif
3330  //Ptcl::vel_cm = stat.pcm.vel;
3331  }
3332 
3333 #ifdef GALPY
3334  std::string galpy_conf_filename = input_parameters.fname_inp.value+".galpy";
3335  galpy_manager.initial(galpy_parameters, stat.time, galpy_conf_filename, restart_flag, print_flag);
3336 #endif
3337 
3338  // set system hard paramters
3339  hard_manager.setDtRange(input_parameters.dt_soft.value/input_parameters.dt_limit_hard_factor.value, input_parameters.dt_min_hermite_index.value);
3340  hard_manager.setEpsSq(input_parameters.eps.value*input_parameters.eps.value);
3341  hard_manager.setGravitationalConstant(input_parameters.gravitational_constant.value);
3342  hard_manager.r_in_base = r_in;
3343  hard_manager.r_out_base = r_out;
3344 #ifdef HARD_CHECK_ENERGY
3345  hard_manager.energy_error_max = input_parameters.e_err_hard.value;
3346 #else
3347  hard_manager.energy_error_max = PS::LARGE_FLOAT;
3348 #endif
3349  hard_manager.n_step_per_orbit = input_parameters.n_step_per_orbit.value;
3350  hard_manager.ap_manager.r_tidal_tensor = r_bin;
3351  hard_manager.ap_manager.id_offset = id_offset;
3352 #ifdef ORBIT_SAMPLING
3353  hard_manager.ap_manager.orbit_manager.setParticleSplitN(input_parameters.n_split.value);
3354 #endif
3355  hard_manager.h4_manager.step.eta_4th = input_parameters.eta.value;
3356  hard_manager.h4_manager.step.eta_2nd = 0.01*input_parameters.eta.value;
3357  hard_manager.h4_manager.step.calcAcc0OffsetSq(mass_average, r_out, input_parameters.gravitational_constant.value);
3358  hard_manager.ar_manager.energy_error_relative_max = input_parameters.e_err_ar.value;
3359  hard_manager.ar_manager.step_count_max = input_parameters.step_limit_ar.value;
3360  hard_manager.ar_manager.step.initialSymplecticCofficients(-6);
3361  hard_manager.ar_manager.slowdown_pert_ratio_ref = input_parameters.sd_factor.value;
3362  hard_manager.ar_manager.slowdown_timescale_max = dt_soft*input_parameters.n_step_per_orbit.value;
3363  //hard_manager.ar_manager.slowdown_timescale_max = dt_soft;
3364 #ifdef SLOWDOWN_MASSRATIO
3365  hard_manager.ar_manager.slowdown_mass_ref = mass_average;
3366 #endif
3367  hard_manager.ar_manager.interrupt_detection_option = input_parameters.interrupt_detection_option.value;
3368 #ifdef STELLAR_EVOLUTION
3369  hard_manager.ar_manager.interaction.stellar_evolution_option = input_parameters.stellar_evolution_option.value;
3370  if (write_style) hard_manager.ar_manager.interaction.stellar_evolution_write_flag = true;
3371  else hard_manager.ar_manager.interaction.stellar_evolution_write_flag = false;
3372 #ifdef BSE_BASE
3373  if (input_parameters.stellar_evolution_option.value>0) {
3374  hard_manager.ar_manager.interaction.bse_manager.initial(bse_parameters, print_flag);
3375  hard_manager.ar_manager.interaction.tide.speed_of_light = hard_manager.ar_manager.interaction.bse_manager.getSpeedOfLight();
3376  }
3377 #endif
3378 #endif
3379 #ifdef ADJUST_GROUP_PRINT
3380  // group information
3381  if (write_style&&input_parameters.adjust_group_write_option.value==1)
3382  hard_manager.h4_manager.adjust_group_write_flag=true;
3383  else
3384  hard_manager.h4_manager.adjust_group_write_flag=false;
3385 #endif
3386 
3387  // check consistence of paramters
3388  input_parameters.checkParams();
3389  hard_manager.checkParams();
3390 
3391  // initial hard class and parameters
3392  system_hard_one_cluster.manager = &hard_manager;
3393  system_hard_one_cluster.setTimeOrigin(stat.time);
3394 
3395  system_hard_isolated.allocateHardIntegrator(input_parameters.n_interrupt_limit.value);
3396  system_hard_isolated.manager = &hard_manager;
3397  system_hard_isolated.setTimeOrigin(stat.time);
3398 
3399 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
3400  system_hard_connected.allocateHardIntegrator(input_parameters.n_interrupt_limit.value);
3401  system_hard_connected.manager = &hard_manager;
3402  system_hard_connected.setTimeOrigin(stat.time);
3403 #endif
3404 
3405  time_kick = stat.time;
3406 
3407  if(write_style>0&&my_rank==0) {
3408  if (print_flag) std::cout<<"----- Dump parameter files -----"<<std::endl;
3409  // save initial parameters
3410  std::string& fname_par = input_parameters.fname_par.value;
3411  if (print_flag) std::cout<<"Save input parameters to file "<<fname_par<<std::endl;
3412  FILE* fpar_out;
3413  if( (fpar_out = fopen(fname_par.c_str(),"w")) == NULL) {
3414  fprintf(stderr,"Error: Cannot open file %s.\n", fname_par.c_str());
3415  abort();
3416  }
3417  input_parameters.input_par_store.writeAscii(fpar_out);
3418  fclose(fpar_out);
3419 
3420  // save hard paramters
3421  std::string fhard_par = input_parameters.fname_par.value + ".hard";
3422  if (print_flag) std::cout<<"Save hard_manager parameters to file "<<fhard_par<<std::endl;
3423  if( (fpar_out = fopen(fhard_par.c_str(),"w")) == NULL) {
3424  fprintf(stderr,"Error: Cannot open file %s.\n", fhard_par.c_str());
3425  abort();
3426  }
3427  hard_manager.writeBinary(fpar_out);
3428  fclose(fpar_out);
3429 
3430 #ifdef BSE_BASE
3431  // save bse parameters
3432  std::string fbse_par = input_parameters.fname_par.value + fbse_par_suffix;
3433  if (print_flag) std::cout<<"Save bse_parameters to file "<<fbse_par<<std::endl;
3434  if( (fpar_out = fopen(fbse_par.c_str(),"w")) == NULL) {
3435  fprintf(stderr,"Error: Cannot open file %s.\n", fbse_par.c_str());
3436  abort();
3437  }
3438  bse_parameters.input_par_store.writeAscii(fpar_out);
3439  fclose(fpar_out);
3440 #endif
3441 
3442 #ifdef GALPY
3443  // save galpy parameters
3444  std::string fgalpy_par = input_parameters.fname_par.value + ".galpy";
3445  if (print_flag) std::cout<<"Save galpy_parameters to file "<<fgalpy_par<<std::endl;
3446  if( (fpar_out = fopen(fgalpy_par.c_str(),"w")) == NULL) {
3447  fprintf(stderr,"Error: Cannot open file %s.\n", fgalpy_par.c_str());
3448  abort();
3449  }
3450  galpy_parameters.input_par_store.writeAscii(fpar_out);
3451  fclose(fpar_out);
3452 #endif
3453  }
3454 
3455  // initial tree step manager
3456  dt_manager.setKDKMode();
3457 
3458  if (print_flag) std::cout<<"----- Finish parameter initialization -----"<<std::endl;
3459 
3460  initial_parameters_flag = true;
3461  }
3462 
3464  void initialStep() {
3465  assert(initial_parameters_flag);
3466  if (initial_step_flag) return;
3467 
3468  assert(checkTimeConsistence());
3469 
3470  // one particle case
3471  if (stat.n_real_glb==1) {
3472  Ptcl::group_data_mode = GroupDataMode::artificial;
3473 
3474  if (stat.n_real_loc==1) system_soft[0].clearForce();
3475 
3477  externalForce();
3478 
3479  // initial status and energy
3480  updateStatus(true);
3481 
3482  PS::F64 dt_tree = input_parameters.dt_soft.value;
3483  dt_manager.setStep(dt_tree);
3484 
3485  // output initial data
3486  file_header.nfile--; // avoid repeating files
3487  output();
3488 
3489 #ifdef PROFILE
3490  clearProfile();
3491 #endif
3492  initial_step_flag = true;
3493 
3494  return;
3495  }
3496 
3497  // domain decomposition
3498  domainDecompose();
3499 
3500  // exchange particles
3501  exchangeParticle();
3502 
3503  PS::F64 dt_tree = input_parameters.dt_soft.value;
3504  dt_manager.setStep(dt_tree);
3505 
3506  // >1. Tree for neighbor searching
3508  treeNeighborSearch();
3509 
3510  // >2. search clusters
3512  Ptcl::group_data_mode = GroupDataMode::cm;
3513  searchCluster();
3514 
3515  // >3. find group and create artificial particles
3517  createGroup(dt_tree);
3518 
3519  // >4 tree soft force
3521  treeSoftForce();
3522 
3524  externalForce();
3525 
3526  // >5 correct change over
3528  treeForceCorrectChangeover();
3529 
3530  // correct force due to the change over update
3531  correctForceChangeOverUpdate();
3532 
3533  // recover mass
3534  kickClusterAndRecoverGroupMemberMass(system_soft, system_hard_isolated.getPtcl(), 0);
3535 
3536 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
3537  // connected
3538  kickClusterAndRecoverGroupMemberMass(system_soft, system_hard_connected.getPtcl(), 0);
3539  // sending list for connected clusters
3540  kickSend(system_soft, search_cluster.getAdrSysConnectClusterSend(), 0);
3541  // send kicked particle from sending list, and receive remote single particle
3542  search_cluster.SendSinglePtcl(system_soft, system_hard_connected.getPtcl());
3543 #endif
3544 
3545  // update global particle system due to kick
3546  writeBackHardParticles();
3547 
3548  // initial status and energy
3549  updateStatus(true);
3550 
3551  // output initial data
3552  file_header.nfile--; // avoid repeating files
3553  output();
3554 
3555  // remove artificial particles
3556  system_soft.setNumberOfParticleLocal(stat.n_real_loc);
3557 
3558 #ifdef CLUSTER_VELOCITY
3559  setParticleGroupDataToCMData();
3560 #endif
3561 
3562 #ifdef PROFILE
3563  clearProfile();
3564 #endif
3565  initial_step_flag = true;
3566  }
3567 
3569  PS::S32 integrateOneToTime(const PS::F64 _time_break=0.0) {
3570  // ensure it is initialized
3571  assert(initial_step_flag);
3572  assert(stat.n_real_glb==1);
3573 
3574  // check time break
3575  PS::F64 time_break = _time_break==0.0? input_parameters.time_end.value: std::min(_time_break,input_parameters.time_end.value);
3576  if (stat.time>=time_break) return 0;
3577 
3578  //PS::F64 dt = std::min(input_parameters.dt_soft.value, time_break-stat.time);
3579  PS::F64 dt_output = input_parameters.dt_snap.value;
3580  //bool start_flag=true;
3581 
3582  if (stat.n_real_loc==1) {
3583  assert(system_soft.getNumberOfParticleLocal()==1);
3584  auto& p = system_soft[0];
3585 
3586  while (true) {
3587 #ifdef PROFILE
3588  profile.total.start();
3589 #endif
3590 
3591  // reset total potential
3592  p.clearForce();
3593 
3595  externalForce();
3596 
3597 #ifdef RECORD_CM_IN_HEADER
3598  stat.calcAndShiftCenterOfMass(&p, stat.n_real_loc);
3599 #endif
3600 
3601  bool interrupt_flag = false; // for interrupt integration when time reach end
3602  bool output_flag = false; // for output snapshot and information
3603 
3604  PS::F64 dt_kick, dt_drift;
3605  // for initial the system
3606  if (dt_manager.isNextStart()) {
3607 
3608  // set step to the begining step
3609  dt_kick = dt_manager.getDtStartContinue();
3610  }
3611  else {
3612  // increase loop counter
3613  n_loop++;
3614 
3615  // for next kick-drift pair
3616  dt_manager.nextContinue();
3617 
3618  if (dt_manager.isNextEndPossible()) {
3619 
3620  // output step, get last kick step
3621  output_flag = (fmod(stat.time, dt_output) == 0.0);
3622 
3623  // check interruption
3624  interrupt_flag = (stat.time>=time_break);
3625 
3626  if (output_flag||interrupt_flag) dt_kick = dt_manager.getDtEndContinue();
3627  else dt_kick = dt_manager.getDtKickContinue();
3628  }
3629  else dt_kick = dt_manager.getDtKickContinue();
3630  }
3631 
3632  //kick
3633  p.vel += p.acc * dt_kick;
3634 #ifdef GALPY
3635  galpy_manager.kickMovePot(dt_kick);
3636 #endif
3637  time_kick += dt_kick;
3638 
3639  // output information
3640  if(output_flag) {
3641 
3642  updateStatus(false);
3643  output();
3644 #ifdef PROFILE
3645  // profile
3646  profile.total.barrier();
3647  PS::Comm::barrier();
3648  profile.total.end();
3649 
3650  printProfile();
3651  clearProfile();
3652 
3653  PS::Comm::barrier();
3654  profile.total.start();
3655 #endif
3656  }
3657 
3658  // interrupt
3659  if(interrupt_flag) {
3660  assert(checkTimeConsistence());
3661 #ifdef PROFILE
3662  profile.total.barrier();
3663  PS::Comm::barrier();
3664  profile.total.end();
3665 #endif
3666  return 0;
3667  }
3668 
3669  // second kick if output exists
3670  if(output_flag) {
3671  dt_kick = dt_manager.getDtStartContinue();
3672  p.vel += p.acc * dt_kick;
3673 #ifdef GALPY
3674  galpy_manager.kickMovePot(dt_kick);
3675 #endif
3676 #ifdef RECORD_CM_IN_HEADER
3677  // correct Ptcl:vel_cm
3678  correctPtclVelCM(dt_kick);
3679 #endif
3680  time_kick += dt_kick;
3681  }
3682 
3683  // get drift step
3684  dt_drift = dt_manager.getDtDriftContinue();
3685 
3686  p.pos += p.vel * dt_drift;
3687 
3688  // drift cm
3689  stat.pcm.pos += stat.pcm.vel*dt_drift;
3690 
3691 #ifdef GALPY
3692  galpy_manager.driftMovePot(dt_drift);
3693 #endif
3694 
3695 #ifdef STELLAR_EVOLUTION
3696  PS::F64 mbk = p.mass;
3697  PS::F64vec vbk = p.vel; //back up velocity in case of change
3698  int modify_flag = hard_manager.ar_manager.interaction.modifyOneParticle(p, stat.time, stat.time + dt_drift);
3699  if (modify_flag) {
3700  auto& v = p.vel;
3701  PS::F64 de_kin = 0.5*(p.mass*(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]) - mbk*(vbk[0]*vbk[0]+vbk[1]*vbk[1]+vbk[2]*vbk[2]));
3702  stat.energy.de_sd_change_cum += de_kin;
3703  stat.energy.de_sd_change_modify_single += de_kin;
3704  stat.energy.de_change_cum += de_kin;
3705  stat.energy.de_change_modify_single += de_kin;
3706  stat.energy.etot_ref += de_kin;
3707  stat.energy.etot_sd_ref += de_kin;
3708  }
3709  // shift time interrupt in order to get consistent time for stellar evolution in the next drift
3710  //p.time_record -= dt;
3711  //p.time_interrupt -= dt;
3712 #endif
3713  stat.time += dt_drift;
3714 
3715 #ifdef PROFILE
3716  profile.total.barrier();
3717  PS::Comm::barrier();
3718  profile.total.end();
3719 
3720  calcProfile();
3721 #endif
3722  }
3723  }
3724  else {
3725  PS::F64 dt_tree = dt_manager.getStep();
3726  while (stat.time < time_break) {
3727  stat.time += dt_tree;
3728  time_kick = stat.time;
3729  }
3730  }
3731 
3732  return 0;
3733  }
3734 
3735 
3737 
3740  PS::S32 integrateToTime(const PS::F64 _time_break=0.0) {
3741 
3742  // ensure it is initialized
3743  assert(initial_step_flag);
3744 
3745  // for one particle case
3746  if (stat.n_real_glb==1) return integrateOneToTime(_time_break);
3747 
3748  // finish interrupted integrations
3749  if (n_interrupt_glb>0) finishInterruptDrift();
3750  // if interrupt still exist, do not continue
3751  if (n_interrupt_glb>0) return n_interrupt_glb;
3752 
3753  // check time break
3754  PS::F64 time_break = _time_break==0.0? input_parameters.time_end.value: std::min(_time_break,input_parameters.time_end.value);
3755  if (stat.time>=time_break) return 0;
3756 
3757  PS::F64 dt_output = input_parameters.dt_snap.value;
3758  PS::F64 dt_tree = dt_manager.getStep();
3759 
3761  while(true) {
3762 #ifdef PROFILE
3763  profile.total.start();
3764 #endif
3765 
3766 #ifdef STELLAR_EVOLUTION
3767  // correct soft potential energy due to mass change
3768  correctSoftPotMassChange();
3769 #endif
3770 
3771  // remove artificial and ununsed particles in system_soft.
3772  removeParticles();
3773 
3774 #ifdef RECORD_CM_IN_HEADER
3775  // update center
3776  stat.calcAndShiftCenterOfMass(&system_soft[0], stat.n_real_loc);
3777 #endif
3778 
3779  // >9. Domain decomposition
3780  domainDecompose();
3781 
3782  // >10. exchange particles
3783  exchangeParticle();
3784 
3785  // >1. Tree for neighbor searching
3787  treeNeighborSearch();
3788 
3789  // >2. search clusters
3791  searchCluster();
3792 
3793  // >3. find group and create artificial particles
3795  createGroup(dt_tree);
3796 
3797  // >4 tree soft force
3799  treeSoftForce() ;
3800 
3802  externalForce();
3803 
3804  // >5 correct change over and potential energy due to mass change
3807  treeForceCorrectChangeover();
3808 
3809 
3810 #ifdef KDKDK_4TH
3811  // only do correction at middle step
3812  if (dt_manager.getCountContinue() == 1) GradientKick();
3813 #endif
3814 
3815  bool interrupt_flag = false; // for interrupt integration when time reach end
3816  bool output_flag = false; // for output snapshot and information
3817  //bool dt_mod_flag = false; // for check whether tree time step need update
3818  bool changeover_flag = false; // for check whether changeover need update
3819  PS::F64 dt_kick, dt_drift;
3820 
3821  // for initial the system
3822  if (dt_manager.isNextStart()) {
3823 
3824  // update changeover if last time it is modified.
3825  correctForceChangeOverUpdate();
3826 
3827  // set step to the begining step
3828  dt_kick = dt_manager.getDtStartContinue();
3829  }
3830  else {
3831 #ifdef PROFILE
3832  profile.other.start();
3833 #endif
3834  // increase loop counter
3835  n_loop++;
3836 
3837  // for next kick-drift pair
3838  dt_manager.nextContinue();
3839 
3840  // update stat time
3841  stat.time = system_hard_one_cluster.getTimeOrigin();
3842 
3843  // check whether output or changeover change are needed (only at the ending step)
3844  if (dt_manager.isNextEndPossible()) {
3845 
3846  // adjust tree step
3847  //dt_mod_flag = adjustDtTreeReduce(dt_reduce_factor, dt_tree, dt_manager.getStep());
3848  //if(dt_mod_flag) dt_kick = dt_manager.getDtEndContinue();
3849 
3850  // output step, get last kick step
3851  output_flag = (fmod(stat.time, dt_output) == 0.0);
3852 
3853  // check changeover change
3854  changeover_flag = (system_hard_isolated.getNClusterChangeOverUpdate()>0);
3855 
3856 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
3857  PS::S32 n_changeover_modify_local = system_hard_connected.getNClusterChangeOverUpdate() + system_hard_isolated.getNClusterChangeOverUpdate();
3858  PS::S32 n_changeover_modify_global = PS::Comm::getSum(n_changeover_modify_local);
3859  if (n_changeover_modify_global>0) changeover_flag = true;
3860 #endif
3861 
3862  // check interruption
3863  interrupt_flag = (stat.time>=time_break);
3864 
3865  // set next step to be last
3866  if (output_flag||changeover_flag||interrupt_flag) dt_kick = dt_manager.getDtEndContinue();
3867  else dt_kick = dt_manager.getDtKickContinue();
3868  }
3869  else dt_kick = dt_manager.getDtKickContinue();
3870 #ifdef PROFILE
3871  profile.other.barrier();
3872  PS::Comm::barrier();
3873  profile.other.end();
3874 #endif
3875  }
3876 
3877 
3878  // >6. kick
3879  kick(dt_kick);
3880 
3881  // >7. write back data
3882  if(output_flag||interrupt_flag) {
3883  // update global particle system due to kick
3884  writeBackHardParticles();
3885  }
3886 
3887  // output information
3888  if(output_flag) {
3889  // update status
3890  updateStatus(false);
3891  output();
3892 
3893 #ifdef PROFILE
3894  // profile
3895  profile.total.barrier();
3896  PS::Comm::barrier();
3897  profile.total.end();
3898 
3899  printProfile();
3900  clearProfile();
3901 
3902  PS::Comm::barrier();
3903  profile.total.start();
3904 #endif
3905  }
3906 
3907  // modify the tree step
3908  //if(dt_mod_flag) {
3909  // dt_manager.setStep(dt_tree/dt_reduce_factor);
3910  // if (input_parameters.print_flag) {
3911  // std::cout<<"Tree time step change, time = "<<stat.time
3912  // <<" dt_soft = "<<dt_tree
3913  // <<" reduce factor = "<<dt_reduce_factor
3914  // <<std::endl;
3915  // }
3916  //}
3917 
3918  // interrupt
3919  if(interrupt_flag) {
3920 #ifdef CLUSTER_VELOCITY
3921  setParticleGroupDataToCMData();
3922 #endif
3923  // correct force due to the change over update
3924  correctForceChangeOverUpdate();
3925 
3926  // need remove artificial particles
3927  system_soft.setNumberOfParticleLocal(stat.n_real_loc);
3928 
3929  assert(checkTimeConsistence());
3930 
3931 #ifdef PROFILE
3932  profile.total.barrier();
3933  PS::Comm::barrier();
3934  profile.total.end();
3935 #endif
3936  return 0;
3937  }
3938 
3939  // second kick if output exists or changeover is modified
3940  //if(dt_mod_flag||output_flag||changeover_flag) {
3941  if(output_flag||changeover_flag) {
3942 
3943  assert(checkTimeConsistence());
3944 
3945  // determine the next tree time step
3946  //PS::F64 dt_reduce_factor_org = 1.0;
3947  //dt_reduce_factor_org = PS::Comm::getMaxValue(dt_reduce_factor_org);
3948  //dt_reduce_factor = 1.0;
3949  //while(dt_reduce_factor<dt_reduce_factor_org) dt_reduce_factor *=2.0;
3950 
3951  //update new tree step if reduce factor is changed
3952  dt_kick = dt_manager.getDtStartContinue();
3953 
3954  correctForceChangeOverUpdate();
3955 
3956  kick(dt_kick);
3957 
3958  }
3959 
3960 
3961  // >8. Hard integration
3962  // get drift step
3963  dt_drift = dt_manager.getDtDriftContinue();
3964 
3965 #ifdef STELLAR_EVOLUTION
3966 #ifdef BSE_BASE
3967  hard_manager.ar_manager.interaction.time_interrupt_max = stat.time + dt_drift;
3968 #endif
3969 #endif
3970 
3971  drift(dt_drift);
3972 
3973  // update stat time
3974  stat.time = system_hard_one_cluster.getTimeOrigin();
3975 
3976 #ifdef PROFILE
3977  // calculate profile
3978  profile.total.barrier();
3979  PS::Comm::barrier();
3980  profile.total.end();
3981 
3982  calcProfile();
3983 #endif
3984 
3985  // when interrupt exist, quit the loop
3986  if (n_interrupt_glb>0) {
3987 #ifdef HARD_INTERRUPT_PRINT
3988  std::cerr<<"Interrupt detected, number: "<<n_interrupt_glb<<std::endl;
3989 #endif
3990  return n_interrupt_glb;
3991  }
3992 
3993  }
3994 
3995  return 0;
3996  }
3997 
3998  void clear() {
3999 
4000  if (fstatus.is_open()) fstatus.close();
4001  if (fesc.is_open()) fesc.close();
4002 #ifdef PROFILE
4003  if (fprofile.is_open()) fprofile.close();
4004 #endif
4005 
4006 #ifdef BSE_BASE
4007  auto& interaction = hard_manager.ar_manager.interaction;
4008  if (interaction.fout_sse.is_open()) interaction.fout_sse.close();
4009  if (interaction.fout_bse.is_open()) interaction.fout_bse.close();
4010 #endif
4011 #ifdef ADJUST_GROUP_PRINT
4012  if (hard_manager.h4_manager.fgroup.is_open()) hard_manager.h4_manager.fgroup.close();
4013 #endif
4014  if (pos_domain) {
4015  delete[] pos_domain;
4016  pos_domain=NULL;
4017  }
4018 
4019  //if (initial_fdps_flag) PS::Finalize();
4020  remove_list.resizeNoInitialize(0);
4021  n_interrupt_glb = 0;
4022  //initial_fdps_flag = false;
4023  read_parameters_flag = false;
4024  read_data_flag = false;
4025  initial_parameters_flag = false;
4026  initial_step_flag = false;
4027  }
4028 
4029  ~PeTar() { clear();}
4030 };
4031 
4032 bool PeTar::initial_fdps_flag = false;
ArtificialParticleManager::getCMParticles
Tptcl * getCMParticles(Tptcl *_ptcl_list)
get c.m. particle list address from a artificial particle array
Definition: artificial_particles.hpp:404
SysCounts::clear
void clear()
Definition: profile.hpp:533
SystemHard::finishIntegrateInterruptClustersOMP
PS::S32 finishIntegrateInterruptClustersOMP()
Finish interrupt integration.
Definition: hard.hpp:2425
EnergyAndMomemtum::ekin
PS::F64 ekin
Definition: energy.hpp:7
Status::ParticleCM::vel
PS::F64vec vel
Definition: status.hpp:21
KickDriftStep::getDtEndContinue
PS::F64 getDtEndContinue()
Get ending step size for continue case.
Definition: kickdriftstep.hpp:161
FileHeader::n_body
long long int n_body
Definition: io.hpp:213
Status::pcm
struct Status::ParticleCM pcm
IOParamsPeTar::dt_soft
IOParams< PS::F64 > dt_soft
Definition: petar.hpp:110
EnergyAndMomemtum::clear
void clear()
Definition: energy.hpp:36
FPSoft
Definition: soft_ptcl.hpp:26
SysProfile::tree_nb
Tprofile tree_nb
Definition: profile.hpp:320
hard.hpp
SystemHard::correctForceWithCutoffClusterOMP
void correctForceWithCutoffClusterOMP(Tsys &_sys, const bool _acorr_flag=false)
Soft force correction due to different cut-off function.
Definition: hard.hpp:3181
Escaper::isEscaper
bool isEscaper(Tptcl &_p, Tpcm &_pcm)
check escaper based on distance and c.m. of the system
Definition: escaper.hpp:14
HardManager::setEpsSq
void setEpsSq(const PS::F64 _eps_sq)
set softening
Definition: hard.hpp:41
HardEnergy::de_sd_change_binary_interrupt
PS::F64 de_sd_change_binary_interrupt
Definition: hard.hpp:126
IOParamsPeTar::r_out
IOParams< PS::F64 > r_out
Definition: petar.hpp:126
IOParamsPeTar::checkParams
bool checkParams()
check paramters
Definition: petar.hpp:583
HardManager::ap_manager
ArtificialParticleManager ap_manager
Definition: hard.hpp:33
FDPSProfile::dump
void dump(std::ostream &fout, const PS::S64 n_loop=1, const PS::S32 width=PROFILE_PRINT_WIDTH)
Definition: profile.hpp:293
SystemHard::updateTimeWriteBack
void updateTimeWriteBack()
Definition: hard.hpp:1969
SystemSoft
PS::ParticleSystem< FPSoft > SystemSoft
Definition: format_transfer.cxx:8
GalpyManager::printReference
static void printReference(std::ostream &fout, const int offset=4)
print reference to cite
Definition: galpy_interface.h:1486
IOParams::value
Type value
Definition: io.hpp:43
SysProfile::output
Tprofile output
Definition: profile.hpp:328
HardEnergy::de_sd
PS::F64 de_sd
Definition: hard.hpp:124
GalpyManager::writePotentialPars
void writePotentialPars(const std::string &_filename, const double &_system_time)
write potential parameters for restart
Definition: galpy_interface.h:1474
BSEManager::printReference
static void printReference(std::ostream &fout, const int offset=4)
print reference to cite
Definition: bse_interface.h:1146
BSEManager::getBSEOutputFilenameSuffix
static std::string getBSEOutputFilenameSuffix()
Definition: bse_interface.h:1178
kickdriftstep.hpp
HardEnergy::de
PS::F64 de
Definition: hard.hpp:120
IOParamsBSE
IO parameters manager for BSE based code.
Definition: bse_interface.h:619
SysCounts::getherClusterCount
void getherClusterCount(int *n, int *count, const long unsigned int size)
Definition: profile.hpp:464
SysCounts::addClusterCount
void addClusterCount(SysCounts &n_count)
Definition: profile.hpp:480
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
domain.hpp
Status::n_escape_glb
PS::S32 n_escape_glb
Definition: status.hpp:15
SysProfile::search_cluster
Tprofile search_cluster
Definition: profile.hpp:324
SysCounts::cluster_isolated
NumCounter cluster_isolated
Definition: profile.hpp:424
IOParamsPeTar::n_leaf_limit
IOParams< PS::S64 > n_leaf_limit
Definition: petar.hpp:95
SearchCluster::n_ptcl_in_multi_cluster_isolated_
PS::ReallocatableArray< PS::S32 > n_ptcl_in_multi_cluster_isolated_
Definition: cluster_list.hpp:178
EnergyAndMomemtum::saveEnergyError
void saveEnergyError()
save current energy error
Definition: energy.hpp:301
ParticleDistributionGenerator::makePlummerModel
static void makePlummerModel(const double mass_glb, const long long int n_glb, const long long int n_loc, double *&mass, PS::F64vec *&pos, PS::F64vec *&vel, const double eng=-0.25, const int seed=0)
Definition: particle_distribution_generator.hpp:173
soft_ptcl.hpp
SysCounts::clusterCount
void clusterCount(const PS::S32 n, const PS::S32 ntimes=1)
Definition: profile.hpp:455
Ptcl::group_data_mode
static GroupDataMode group_data_mode
Definition: ptcl.hpp:47
Ptcl::mean_mass_inv
static PS::F64 mean_mass_inv
Definition: ptcl.hpp:45
SysCounts::hard_connected
NumCounter hard_connected
Definition: profile.hpp:422
IOParamsPeTar::r_bin
IOParams< PS::F64 > r_bin
Definition: petar.hpp:127
FileHeader::time
double time
Definition: io.hpp:214
GalpyManager::updatePotential
void updatePotential(const double &_system_time, const bool _print_flag)
Update Potential types and arguments based on time.
Definition: galpy_interface.h:1443
PIKG::F64vec
Vector3< F64 > F64vec
Definition: pikg_vector.hpp:167
SearchCluster::mediator_sorted_id_cluster_
PS::ReallocatableArray< Mediator > mediator_sorted_id_cluster_
Definition: cluster_list.hpp:180
IOParamsContainer::readAscii
void readAscii(FILE *_fin)
Definition: io.hpp:112
IOParamsBSE::input_par_store
IOParamsContainer input_par_store
Definition: bse_interface.h:621
KickDriftStep
Definition: kickdriftstep.hpp:1
SystemHard::getNClusterChangeOverUpdate
PS::S32 getNClusterChangeOverUpdate() const
Definition: hard.hpp:1961
CalcForceEpSpQuadNoSimd
Definition: soft_force.hpp:160
HardManager::r_out_base
PS::F64 r_out_base
Definition: hard.hpp:31
SystemHard::getNumberOfInterruptClusters
PS::S32 getNumberOfInterruptClusters() const
Definition: hard.hpp:1933
IOParamsContainer::printHelp
void printHelp(std::ostream &os, const int _offset_short_key=2, const int _offset_long_key=1, const int _width_key=23) const
Definition: io.hpp:198
IOParamsPeTar::unit_set
IOParams< PS::S64 > unit_set
Definition: petar.hpp:107
galpy_pot_movie.dt
dt
Definition: galpy_pot_movie.py:134
GroupDataMode::cm
@ cm
HardManager::energy_error_max
PS::F64 energy_error_max
Definition: hard.hpp:28
Tprofile::end
void end()
Definition: profile.hpp:191
SysCounts::ARC_n_groups_iso
NumCounter ARC_n_groups_iso
Definition: profile.hpp:429
SysCounts::ep_sp_interact
NumCounter ep_sp_interact
Definition: profile.hpp:433
GalpyManager::initial
void initial(const IOParamsGalpy &_input, const double _time, const std::string _conf_name, const bool _restart_flag, const bool _print_flag=false)
initialization function
Definition: galpy_interface.h:852
Tprofile::barrier
void barrier()
Definition: profile.hpp:187
SysProfile::dumpBarrier
void dumpBarrier(std::ostream &fout, const PS::S64 n_loop=1, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:375
IOParamsPeTar::theta
IOParams< PS::F64 > theta
Definition: petar.hpp:94
HardManager
Hard integrator parameter manager.
Definition: hard.hpp:26
CalcForceEpEpWithLinearCutoffNoSimd
FORCE FUNCTOR.
Definition: soft_force.hpp:38
SearchCluster::initialize
void initialize()
Definition: cluster_list.hpp:333
RetrieveForceCUDA
PS::S32 RetrieveForceCUDA(const PS::S32 tag, const PS::S32 n_walk, const PS::S32 *ni, ForceSoft **force)
IOParamsPeTar::r_search_min
IOParams< PS::F64 > r_search_min
Definition: petar.hpp:129
IOParamsPeTar::n_step_per_orbit
IOParams< PS::S64 > n_step_per_orbit
Definition: petar.hpp:103
FPSoft::pot_tot
PS::F64 pot_tot
Definition: soft_ptcl.hpp:32
IOParamsBSE::vscale
IOParams< double > vscale
Definition: bse_interface.h:663
KickDriftStep::getStep
PS::F64 getStep() const
get base step size
Definition: kickdriftstep.hpp:107
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
SysProfile
Definition: profile.hpp:313
IOParamsPeTar::dt_min_hermite_index
IOParams< PS::S64 > dt_min_hermite_index
Definition: petar.hpp:116
GalpyManager::calcAccPot
void calcAccPot(double *acc, double &pot, const double _time, const double gm, const double *pos_g, const double *pos_l)
calculate acceleration and potential at give position
Definition: galpy_interface.h:1603
HardEnergy::de_change_cum
PS::F64 de_change_cum
Definition: hard.hpp:121
get_version.hpp
IOParamsGalpy::print_flag
bool print_flag
Definition: galpy_interface.h:28
IOParamsPeTar
IO parameters for Petar.
Definition: petar.hpp:89
IOParams::print
void print(std::ostream &os) const
Definition: io.hpp:54
HardManager::h4_manager
H4::HermiteManager< HermiteInteraction > h4_manager
Definition: hard.hpp:34
SysProfile::getMin
SysProfile getMin()
Definition: profile.hpp:399
IOParamsPeTar::IOParamsPeTar
IOParamsPeTar()
Definition: petar.hpp:151
KickDriftStep::setStep
void setStep(const PS::F64 _ds)
set base step size
Definition: kickdriftstep.hpp:76
particle_distribution_generator.hpp
GalpyManager::kickMovePot
void kickMovePot(const double dt)
kick velocity of moving potential (mode 2)
Definition: galpy_interface.h:1503
IOParamsPeTar::data_format
IOParams< PS::S64 > data_format
Definition: petar.hpp:132
SysProfile::hard_isolated
Tprofile hard_isolated
Definition: profile.hpp:317
BSEManager::getBSEName
static std::string getBSEName()
Definition: bse_interface.h:1198
SystemHard::correctPotWithCutoffOMP
void correctPotWithCutoffOMP(Tsys &_sys, const PS::ReallocatableArray< PS::S32 > &_ptcl_list)
potential correction for single cluster
Definition: hard.hpp:3102
force_gpu_cuda.hpp
GroupDataMode::artificial
@ artificial
GalpyManager
A class to manager the API to Galpy.
Definition: galpy_interface.h:658
Status::time
PS::F64 time
Definition: status.hpp:9
SysProfile::force_correct
Tprofile force_correct
Definition: profile.hpp:322
SearchNeighborEpEpNoSimd
Definition: soft_force.hpp:11
Status
class for measure the status of the system
Definition: status.hpp:7
SysProfile::dump
void dump(std::ostream &fout, const PS::S64 n_loop=1, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:361
IOParamsPeTar::append_switcher
IOParams< PS::S64 > append_switcher
Definition: petar.hpp:141
SysProfile::tree_soft
Tprofile tree_soft
Definition: profile.hpp:321
IOParamsGalpy::input_par_store
IOParamsContainer input_par_store
Definition: galpy_interface.h:17
KickDriftStep::isNextStart
bool isNextStart() const
get whether next is start
Definition: kickdriftstep.hpp:172
PIKG::S64
int64_t S64
Definition: pikg_vector.hpp:23
SysCounts::ep_ep_interact
NumCounter ep_ep_interact
Definition: profile.hpp:432
KickDriftStep::isNextEndPossible
bool isNextEndPossible() const
get whether next can be end
Definition: kickdriftstep.hpp:177
SystemHard::setTimeOrigin
void setTimeOrigin(const PS::F64 _time_origin)
Definition: hard.hpp:1965
SysCounts::hard_interrupt
NumCounter hard_interrupt
Definition: profile.hpp:423
SysProfile::kick
Tprofile kick
Definition: profile.hpp:323
EPISoft::r_out
static PS::F64 r_out
Definition: soft_ptcl.hpp:282
IOParamsPeTar::r_escape
IOParams< PS::F64 > r_escape
Definition: petar.hpp:130
ArtificialParticleManager::r_tidal_tensor
PS::F64 r_tidal_tensor
Definition: artificial_particles.hpp:235
IOParams::key
const char * key
Definition: io.hpp:44
IOParamsPeTar::time_end
IOParams< PS::F64 > time_end
Definition: petar.hpp:104
FPSoft::acc
PS::F64vec acc
Definition: soft_ptcl.hpp:28
io.hpp
IOParamsPeTar::write_style
IOParams< PS::S64 > write_style
Definition: petar.hpp:133
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
Ptcl::r_search_min
static PS::F64 r_search_min
Definition: ptcl.hpp:43
ParticleDistributionGenerator::makeKeplerDisk
static void makeKeplerDisk(PS::F64 &mass_planet_glb, PS::F64 *&mass, PS::F64vec *&pos, PS::F64vec *&vel, const long long int n_glb, const long long int n_loc, const double a_in, const double a_out, const double e_rms, const double i_rms, const double dens=10.0, const double mass_sun=1.0, const double a_ice=0.0, const double f_ice=1.0, const double power=-1.5, const int seed=0, const double gravitational_constant=1.0)
Definition: particle_distribution_generator.hpp:111
KickDriftStep::setKDKMode
void setKDKMode()
set DKD mode (default is KDK)
Definition: kickdriftstep.hpp:66
SysProfile::total
Tprofile total
Definition: profile.hpp:315
profile.hpp
SysProfile::domain
Tprofile domain
Definition: profile.hpp:326
Status::n_remove_glb
PS::S32 n_remove_glb
Definition: status.hpp:14
GalpyManager::printData
void printData(std::ostream &fout)
print current potential data
Definition: galpy_interface.h:767
SysCounts::n_cluster
std::map< PS::S32, PS::S32 > n_cluster
Histogram of number of particles in clusters.
Definition: profile.hpp:436
Status::printColumnTitle
void printColumnTitle(std::ofstream &_fout, const PS::S32 _width=20)
print titles of class members using column style
Definition: status.hpp:256
HardEnergy::de_change_binary_interrupt
PS::F64 de_change_binary_interrupt
Definition: hard.hpp:122
HardManager::ar_manager
AR::TimeTransformedSymplecticManager< ARInteraction > ar_manager
Definition: hard.hpp:35
SystemHard::setPtclForOneClusterOMP
void setPtclForOneClusterOMP(const Tsys &sys, const PS::ReallocatableArray< PS::S32 > &adr_array)
Definition: hard.hpp:2026
IOParamsPeTar::search_vel_factor
IOParams< PS::F64 > search_vel_factor
Definition: petar.hpp:113
HardManager::r_in_base
PS::F64 r_in_base
Definition: hard.hpp:30
IOParamsPeTar::gravitational_constant
IOParams< PS::F64 > gravitational_constant
Definition: petar.hpp:106
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
IOParamsPeTar::n_group_limit
IOParams< PS::S64 > n_group_limit
Definition: petar.hpp:96
IOParamsPeTar::e_err_ar
IOParams< PS::F64 > e_err_ar
Definition: petar.hpp:120
IOParamsPeTar::ratio_r_cut
IOParams< PS::F64 > ratio_r_cut
Definition: petar.hpp:93
Status::n_real_glb
PS::S64 n_real_glb
Definition: status.hpp:11
ArtificialParticleManager::orbit_manager
OrbitManager orbit_manager
gravitational constant
Definition: artificial_particles.hpp:238
IOParamsPeTar::n_smp_ave
IOParams< PS::S64 > n_smp_ave
Definition: petar.hpp:98
Status::n_all_glb
PS::S64 n_all_glb
Definition: status.hpp:13
IOParamsPeTar::n_glb
IOParams< PS::S64 > n_glb
Definition: petar.hpp:108
SearchCluster::getAdrSysConnectClusterSend
const PS::ReallocatableArray< PS::S32 > & getAdrSysConnectClusterSend()
get the address list of ptcl need to send in connected clusters
Definition: cluster_list.hpp:1376
IOParamsPeTar::input_par_store
IOParamsContainer input_par_store
Definition: petar.hpp:92
SystemHard::getNumberOfClusters
PS::S32 getNumberOfClusters() const
Definition: hard.hpp:1929
IOParamsPeTar::nstep_dt_soft_kepler
IOParams< PS::F64 > nstep_dt_soft_kepler
Definition: petar.hpp:112
galpy_interface.h
ArtificialParticleManager::getArtificialParticleN
PS::S32 getArtificialParticleN() const
get artificial particle total number
Definition: artificial_particles.hpp:427
Status::n_all_loc
PS::S64 n_all_loc
Definition: status.hpp:12
cluster_list.hpp
HardManager::setGravitationalConstant
void setGravitationalConstant(const PS::F64 _g)
set gravitational constant
Definition: hard.hpp:48
SystemHard
Hard system.
Definition: hard.hpp:1190
escaper.hpp
EnergyAndMomemtum::etot_ref
PS::F64 etot_ref
Definition: energy.hpp:9
IOParamsPeTar::fname_inp
IOParams< std::string > fname_inp
Definition: petar.hpp:144
SysProfile::create_group
Tprofile create_group
Definition: profile.hpp:325
Escaper
escaper checker
Definition: escaper.hpp:5
SearchCluster::adr_sys_multi_cluster_isolated_
PS::ReallocatableArray< PS::S32 > adr_sys_multi_cluster_isolated_
Definition: cluster_list.hpp:182
IOParamsContainer::writeAscii
void writeAscii(FILE *_fout)
Definition: io.hpp:106
Tprofile::tbar
PS::F64 tbar
Definition: profile.hpp:178
KickDriftStep::getCountContinue
PS::S32 getCountContinue() const
get step count for continue case
Definition: kickdriftstep.hpp:51
IOParamsPeTar::dt_snap
IOParams< PS::F64 > dt_snap
Definition: petar.hpp:111
status.hpp
SearchCluster::ptcl_recv_
PS::ReallocatableArray< PtclComm > ptcl_recv_
Definition: cluster_list.hpp:181
Status::print
void print(std::ostream &_fout, const PS::S32 _precision=7)
print title and values in one lines
Definition: status.hpp:290
SysProfile::hard_connected
Tprofile hard_connected
Definition: profile.hpp:318
SysProfile::dumpName
void dumpName(std::ostream &fout, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:368
IOParamsPeTar::n_bin
IOParams< PS::S64 > n_bin
Definition: petar.hpp:102
SysProfile::clear
void clear()
Definition: profile.hpp:409
EPJSoft
Definition: soft_ptcl.hpp:311
Status::ParticleCM::pos
PS::F64vec pos
Definition: status.hpp:20
WRITE_PRECISION
const int WRITE_PRECISION
Definition: bse_test.cxx:11
HardEnergy::epot_sd_correction
PS::F64 epot_sd_correction
Definition: hard.hpp:129
ArtificialParticleManager::id_offset
PS::S64 id_offset
tidal tensor maximum distance of particles
Definition: artificial_particles.hpp:236
CalcForceWithLinearCutoffCUDA
Definition: force_gpu_cuda.hpp:137
SysProfile::hard_interrupt
Tprofile hard_interrupt
Definition: profile.hpp:319
IOParamsBSE::rscale
IOParams< double > rscale
Definition: bse_interface.h:661
FDPSProfile
Definition: profile.hpp:267
FPSoft::adr
PS::S32 adr
Definition: soft_ptcl.hpp:42
SysProfile::other
Tprofile other
Definition: profile.hpp:330
IOParamsPeTar::eps
IOParams< PS::F64 > eps
Definition: petar.hpp:125
IOParamsPeTar::step_limit_ar
IOParams< PS::S64 > step_limit_ar
Definition: petar.hpp:124
SearchCluster::searchClusterLocal
void searchClusterLocal()
Definition: cluster_list.hpp:613
FPSoft::printTitleWithMeaning
static int printTitleWithMeaning(std::ostream &_fout, const int _counter=0, const int _offset=0)
print column title with meaning (each line for one column)
Definition: soft_ptcl.hpp:203
IOParams< PS::F64 >
IOParamsPeTar::sd_factor
IOParams< PS::F64 > sd_factor
Definition: petar.hpp:131
SystemHard::setPtclForIsolatedMultiClusterOMP
void setPtclForIsolatedMultiClusterOMP(const Tsys &sys, const PS::ReallocatableArray< PS::S32 > &_adr_array, const PS::ReallocatableArray< PS::S32 > &_n_ptcl_in_cluster)
Definition: hard.hpp:2229
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
IOParamsPeTar::eta
IOParams< PS::F64 > eta
Definition: petar.hpp:105
SystemHard::writeBackPtclForMultiCluster
void writeBackPtclForMultiCluster(Tsys &_sys, PS::ReallocatableArray< PS::S32 > &_mass_modify_list)
write back hard particles to global system, check mass modification and update time of write back
Definition: hard.hpp:2192
IOParamsGalpy::read
int read(int argc, char *argv[], const int opt_used_pre=0)
reading parameters from GNU option API
Definition: galpy_interface.h:49
SystemHard::allocateHardIntegrator
void allocateHardIntegrator(const PS::S32 _n_hard_int)
allocate memorgy for HardIntegrator
Definition: hard.hpp:1836
SystemHard::getClusterNumberOfMemberList
PS::S32 * getClusterNumberOfMemberList(const std::size_t i=0) const
Definition: hard.hpp:1941
force_fugaku.hpp
EPISoft::eps
static PS::F64 eps
Definition: soft_ptcl.hpp:281
SysCounts::cluster_connected
NumCounter cluster_connected
Definition: profile.hpp:425
WRITE_WIDTH
const int WRITE_WIDTH
Definition: bse_test.cxx:10
IOParamsGalpy
IO parameters for Galpy manager.
Definition: galpy_interface.h:15
SysCounts::hard_single
NumCounter hard_single
Definition: profile.hpp:420
KickDriftStep::getDtKickContinue
PS::F64 getDtKickContinue()
Get kick step size for continue case.
Definition: kickdriftstep.hpp:122
IOParamsPeTar::read
int read(int argc, char *argv[], const int opt_used_pre=0)
reading parameters from GNU option API
Definition: petar.hpp:225
IOParamsPeTar::dt_limit_hard_factor
IOParams< PS::F64 > dt_limit_hard_factor
Definition: petar.hpp:115
SystemHard::correctForceForChangeOverUpdateOMP
void correctForceForChangeOverUpdateOMP(Tsys &_sys, Ttree &_tree, const PS::S32 *_adr_send=NULL, const PS::S32 _n_send=0)
Correct force due to change over factor change.
Definition: hard.hpp:3262
PIKG::max
T max(const Vector3< T > &v)
Definition: pikg_vector.hpp:143
SystemHard::findGroupsAndCreateArtificialParticlesOMP
void findGroupsAndCreateArtificialParticlesOMP(Tsys &_sys, const PS::F64 _dt_tree)
Find groups and create aritfical particles to sys.
Definition: hard.hpp:2757
FileHeader::nfile
long long int nfile
Definition: io.hpp:212
hard_assert.hpp
SysCounts::hard_isolated
NumCounter hard_isolated
Definition: profile.hpp:421
SysProfile::exchange
Tprofile exchange
Definition: profile.hpp:327
HardManager::checkParams
bool checkParams()
check paramters
Definition: hard.hpp:68
SysProfile::dumpBarrierName
void dumpBarrierName(std::ostream &fout, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:382
IOParamsPeTar::fname_par
IOParams< std::string > fname_par
Definition: petar.hpp:143
IOParamsPeTar::n_interrupt_limit
IOParams< PS::S64 > n_interrupt_limit
Definition: petar.hpp:97
static_variables.hpp
SysCounts::clearClusterCount
void clearClusterCount()
Definition: profile.hpp:460
HardEnergy::de_sd_change_modify_single
PS::F64 de_sd_change_modify_single
Definition: hard.hpp:127
KickDriftStep::getDtStartContinue
PS::F64 getDtStartContinue()
Get staring step size for continue case.
Definition: kickdriftstep.hpp:150
Status::n_real_loc
PS::S64 n_real_loc
Definition: status.hpp:10
SearchCluster
Definition: cluster_list.hpp:176
SystemHard::manager
HardManager * manager
Definition: hard.hpp:1221
HardEnergy::de_change_modify_single
PS::F64 de_change_modify_single
Definition: hard.hpp:123
IOParamsPeTar::interrupt_detection_option
IOParams< PS::S64 > interrupt_detection_option
Definition: petar.hpp:137
FDPSProfile::dumpName
void dumpName(std::ostream &fout, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:274
SysCounts::dumpName
void dumpName(std::ostream &fout, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:526
SystemHard::initializeForOneCluster
void initializeForOneCluster(const PS::S32 n)
Definition: hard.hpp:1851
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
Status::energy
EnergyAndMomemtum energy
Definition: status.hpp:17
GalpyManager::driftMovePot
void driftMovePot(const double dt)
drift position of moving potential (mode 2)
Definition: galpy_interface.h:1522
Escaper::check_energy_flag
bool check_energy_flag
Definition: escaper.hpp:8
IOParamsBSE::mscale
IOParams< double > mscale
Definition: bse_interface.h:662
KickDriftStep::getDtDriftContinue
PS::F64 getDtDriftContinue()
Get drift step size for continue case.
Definition: kickdriftstep.hpp:136
ArtificialParticleManager
class to organize artificial particles
Definition: artificial_particles.hpp:232
GalpyManager::resetPotAcc
void resetPotAcc()
reset acceleraltion of potential
Definition: galpy_interface.h:1492
HardManager::setDtRange
void setDtRange(const PS::F64 _dt_max, const PS::S32 _dt_min_index)
set time step range
Definition: hard.hpp:61
Escaper::r_escape_sq
PS::F64 r_escape_sq
Definition: escaper.hpp:7
BSEManager::getSSEOutputFilenameSuffix
static std::string getSSEOutputFilenameSuffix()
Definition: bse_interface.h:1168
SysCounts::ARC_substep_sum
NumCounter ARC_substep_sum
Definition: profile.hpp:426
HardEnergy::de_sd_change_cum
PS::F64 de_sd_change_cum
Definition: hard.hpp:125
Status::ParticleCM::is_center_shift_flag
bool is_center_shift_flag
Definition: status.hpp:22
SystemHard::getTimeOrigin
PS::F64 getTimeOrigin() const
Definition: hard.hpp:1973
SystemHard::getPtcl
PS::ReallocatableArray< PtclH4 > & getPtcl()
Definition: hard.hpp:1925
FileHeader
Definition: io.hpp:210
Ptcl::search_factor
static PS::F64 search_factor
Definition: ptcl.hpp:42
energy.hpp
Tprofile::time
PS::F64 time
Definition: profile.hpp:177
PRINT_PRECISION
#define PRINT_PRECISION
Definition: io.hpp:11
Status::calcCenterOfMass
void calcCenterOfMass(Tsoft *_tsys, const PS::S64 _n, int _mode=3)
calculate the center of system
Definition: status.hpp:116
SystemHard::writeBackPtclForOneClusterOMP
void writeBackPtclForOneClusterOMP(Tsys &sys, PS::ReallocatableArray< PS::S32 > &_mass_modify_list)
write back hard particles to global system and update time of write back OMP version
Definition: hard.hpp:2143
SysCounts::H4_step_sum
NumCounter H4_step_sum
Definition: profile.hpp:430
SystemHard::correctForceWithCutoffTreeNeighborAndClusterOMP
void correctForceWithCutoffTreeNeighborAndClusterOMP(Tsys &_sys, Ttree &_tree, const PS::ReallocatableArray< PS::S32 > &_adr_send, const bool _acorr_flag=false)
Soft force correction due to different cut-off function.
Definition: hard.hpp:3132
SysProfile::hard_single
Tprofile hard_single
Definition: profile.hpp:316
Ptcl::id
PS::S64 id
Definition: ptcl.hpp:39
Ptcl::r_group_crit_ratio
static PS::F64 r_group_crit_ratio
Definition: ptcl.hpp:44
IOParamsContainer
IO Params container.
Definition: io.hpp:82
GetVersion
std::string GetVersion()
Definition: get_version.hpp:4
SysCounts::printHist
void printHist(std::ostream &fout, const PS::S64 n_loop=1, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:487
IOParamsBSE::tscale
IOParams< double > tscale
Definition: bse_interface.h:660
Status::ParticleCM::mass
PS::F64 mass
Definition: status.hpp:19
CalcForceEpSpMonoNoSimd
Definition: soft_force.hpp:125
IOParamsPeTar::search_peri_factor
IOParams< PS::F64 > search_peri_factor
Definition: petar.hpp:114
IOParamsPeTar::id_offset
IOParams< PS::S64 > id_offset
Definition: petar.hpp:109
SystemHard::driveForOneClusterOMP
void driveForOneClusterOMP(const PS::F64 _dt)
integrate one isolated particle
Definition: hard.hpp:2046
PIKG::U64
uint64_t U64
Definition: pikg_vector.hpp:20
SysCounts::dump
void dump(std::ostream &fout, const PS::S64 n_loop=1, const PS::S32 width=PROFILE_PRINT_WIDTH) const
Definition: profile.hpp:513
IOParamsBSE::read
int read(int argc, char *argv[], const int opt_used_pre=0)
reading parameters from GNU option API
Definition: bse_interface.h:754
HardEnergy
Definition: hard.hpp:119
SearchCluster::setIdClusterLocal
void setIdClusterLocal()
Definition: cluster_list.hpp:677
HardManager::n_step_per_orbit
PS::F64 n_step_per_orbit
Definition: hard.hpp:32
Status::printColumn
void printColumn(std::ofstream &_fout, const PS::S32 _width=20)
print data of class members using column style
Definition: status.hpp:273
SearchCluster::getAdrSysOneCluster
const PS::ReallocatableArray< PS::S32 > & getAdrSysOneCluster()
Definition: cluster_list.hpp:1371
SysProfile::status
Tprofile status
Definition: profile.hpp:329
PIKG::F32
float F32
Definition: pikg_vector.hpp:18
SysCounts
Definition: profile.hpp:418
GalpyManager::calcMovePotAccFromPot
void calcMovePotAccFromPot(const double _time, const double *pos_pcm)
calculate acceleration for moving potentials from other potentials
Definition: galpy_interface.h:1543
SearchCluster::searchNeighborOMP
void searchNeighborOMP(Tsys &sys, Ttree &tree, const PS::F64ort pos_domain[], const PS::F64 _G, const PS::F64 _radius_factor)
search neighbors and separate isolated and multiple clusters
Definition: cluster_list.hpp:425
soft_force.hpp
SysCounts::n_neighbor_zero
NumCounter n_neighbor_zero
Definition: profile.hpp:431
ForceSoft::grav_const
static PS::F64 grav_const
neighbor number+1
Definition: soft_ptcl.hpp:15
SystemHard::setPtclForConnectedCluster
void setPtclForConnectedCluster(const Tsys &sys, const PS::ReallocatableArray< Tmediator > &med, const PS::ReallocatableArray< Tptcl > &ptcl_recv)
Definition: hard.hpp:1861
KickDriftStep::nextContinue
void nextContinue()
advance step count for continuing case
Definition: kickdriftstep.hpp:34
HardManager::writeBinary
void writeBinary(FILE *_fp)
write class data to file with binary format
Definition: hard.hpp:83
HardEnergy::ekin_sd_correction
PS::F64 ekin_sd_correction
Definition: hard.hpp:128
IOParamsPeTar::print_flag
bool print_flag
Definition: petar.hpp:147
EnergyAndMomemtum::epot
PS::F64 epot
Definition: energy.hpp:8
IOParamsPeTar::update_changeover_flag
bool update_changeover_flag
Definition: petar.hpp:148
IOParamsPeTar::fname_snp
IOParams< std::string > fname_snp
Definition: petar.hpp:142
SystemHard::driveForMultiClusterOMP
int driveForMultiClusterOMP(const PS::F64 dt, Tpsoft *_ptcl_soft)
Hard integration for clusters.
Definition: hard.hpp:2271
SysCounts::ARC_tsyn_step_sum
NumCounter ARC_tsyn_step_sum
Definition: profile.hpp:427
IOParamsPeTar::update_rsearch_flag
bool update_rsearch_flag
Definition: petar.hpp:149
SysCounts::ARC_n_groups
NumCounter ARC_n_groups
Definition: profile.hpp:428
Tprofile::start
void start()
Definition: profile.hpp:183
IOParamsBSE::print_flag
bool print_flag
Definition: bse_interface.h:666