flat_hermite.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\     
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\
00010 
00011 // hermite.C
00012 // hermite scheme with simple (flat-tree) hdyn structure
00013 //
00014 // J. Makino    92-12-04:  Seems to work... Energy conservation looks good.
00015 // S. McMillan  93-04-02:  Cleaned up and modified command-line interface.
00016 // J. Makino    96-07-xx:  Seems not to work... Energy conservation looks bad.
00017 
00018 void hdyn::flat_set_first_timestep(real eta_for_firststep, real max_step_size)
00019 {
00020     real a1scaled2 = jerk * jerk;
00021     real a2 = acc * acc;
00022     real newstep;
00023 
00024     if (a1scaled2 > 0.0) {
00025         newstep = eta_for_firststep * sqrt(a2 / a1scaled2);
00026     } else {
00027         newstep = eta_for_firststep * 0.125;
00028     }
00029     newstep = adjust_number_to_power(newstep, max_step_size);
00030     timestep = newstep;
00031 }
00032 
00033 real flat_new_timestep(vector & at3,    // third order term
00034                        vector & bt2,    // 2nd order term
00035                        vector & jerk,   // 1st order term
00036                        vector & acc,    // 0th order term
00037                        real timestep,   // old timestep
00038                        real time,       // present time
00039                        real eta,        // accuracy parameter
00040                        real dtmax)      // maximum stepsize
00041 {
00042     real a3scaled2 = 36 * at3 * at3;
00043     real a2scaled2 = 4 * bt2 * bt2;
00044     real a1scaled2 = timestep * jerk * jerk;
00045     real a2 = acc * acc;
00046 
00047 // simple criterion:
00048     real newstep = eta * sqrt(sqrt(a2 / a2scaled2)) * timestep;
00049 
00050     if (newstep < timestep) {
00051         return timestep * 0.5;
00052     } else if (newstep < 2 * timestep) {
00053         return timestep;
00054     } else if (fmod(time, 2 * timestep) != 0.0 || 2 * timestep > dtmax) {
00055         return timestep;
00056     } else {
00057         return timestep * 2;
00058     }
00059 }
00060 
00061 void hdyn::flat_update(const real eta, const real dtmax)
00062 {
00063     vector at3 = 2 * (old_acc - acc) + timestep * (old_jerk + jerk);
00064     vector bt2 = -3 * (old_acc - acc) - timestep * (2 * old_jerk + jerk);
00065     time = time + timestep;
00066     timestep = flat_new_timestep(at3, bt2, jerk, acc,
00067                                  timestep, time, eta, dtmax);
00068 }
00069 
00070 void accumulate_energies(hdyn * b, real & epot, real & ekin, real & etot)
00071 {
00072     if (b->get_oldest_daughter() != NULL) {
00073         for (hdyn * bb = b->get_oldest_daughter(); bb != NULL;
00074              bb = bb->get_younger_sister()) {
00075             accumulate_energies(bb, epot, ekin, etot);
00076         }
00077     } else {
00078         real mi = b->get_mass();
00079         epot += 0.5 * mi * b->get_pot();
00080         vector vel = b->get_vel();
00081         ekin += 0.5 * mi * vel * vel;
00082     }
00083     etot = ekin + epot;
00084 }
00085 
00086 void compute_energies(hdyn * b, real & epot, real & ekin, real & etot)
00087 {
00088     epot = ekin = etot = 0;
00089     accumulate_energies(b, epot, ekin, etot);
00090 }
00091 
00092 void print_energies(hdyn * b, int mode)
00093 {
00094     real epot;
00095     real ekin;
00096     real etot;
00097     static real old_etot;
00098     compute_energies(b, epot, ekin, etot);
00099     cerr << "    Energies:  " << ekin << "  " << epot << "  " << etot;
00100     if (mode) {
00101         real de = (etot - old_etot) / etot;
00102         cerr << ";  de = " << de << "\n";
00103     } else {
00104         old_etot = etot;
00105         cerr << "\n";
00106     }
00107     hdyn *dummy = b;            // to make the compiler happy
00108 
00109 }
00110 
00111 hdyn *particle_to_move(hdyn * b, real & tmin)
00112 {
00113     hdyn *bmin = NULL;
00114     if (b->get_oldest_daughter() != NULL) {
00115         for (hdyn * bb = b->get_oldest_daughter(); bb != NULL;
00116              bb = bb->get_younger_sister()) {
00117             hdyn *btmp;
00118             real ttmp = 1e100;
00119             btmp = particle_to_move(bb, ttmp);
00120             if (ttmp < tmin) {
00121                 tmin = ttmp;
00122                 bmin = btmp;
00123             }
00124         }
00125     } else {
00126         if (tmin > b->get_next_time()) {
00127             bmin = b;
00128             tmin = b->get_next_time();
00129         }
00130     }
00131     return bmin;
00132 }
00133 
00134 void flat_initialize_system_phase2(hdyn * b, hdyn * root, real eps,
00135                                    real eta, real dtout)
00136 {
00137     if (b->get_oldest_daughter() != NULL) {
00138         for (hdyn * bb = b->get_oldest_daughter(); bb != NULL;
00139              bb = bb->get_younger_sister()) {
00140             flat_initialize_system_phase2(bb, root, eps, eta, dtout);
00141         }
00142     } else {
00143         b->clear_interaction();
00144         b->flat_calculate_acc_and_jerk(root, eps * eps);
00145         b->store_old_force();
00146         b->flat_set_first_timestep(0.5 * eta, dtout);
00147     }
00148 }
00149 
00150 local void shift_cm(hdyn * b, real snap_cube_size)
00151 {
00152     real mass = 0;
00153     vector cmpos = vector(0, 0, 0);
00154 
00155     for_all_daughters(hdyn, b, bb) {
00156         vector x = bb->get_pos();
00157         if (abs(x[0]) < snap_cube_size
00158             && abs(x[1]) < snap_cube_size
00159             && abs(x[2]) < snap_cube_size) {
00160             mass += bb->get_mass();
00161             cmpos += bb->get_mass() * bb->get_pos();
00162         }
00163     }
00164 
00165     if (mass > 0) {
00166         cmpos /= mass;
00167         for_all_daughters(hdyn, b, bbb) bbb->inc_pos(-cmpos);
00168     }
00169 }
00170 
00171 void evolve_system(hdyn * b,            // hdyn array                   
00172                    real delta_t,        // time span of the integration 
00173                    real eta,            // accuracy parameter 
00174                    real dtout,          // output time interval
00175                    real dt_snap,        // snapshot output interval
00176                    real eps,            // softening length             
00177                    real snap_cube_size)
00178 {
00179     real t = b->get_time();     // current time
00180     real t_end = t + delta_t;   // final time, at the end of the integration
00181     real t_next = t + dtout;    // time of next output
00182     real t_snap = t + dt_snap;  // time of next snapshot;
00183 
00184     int steps = 0;
00185 
00186     initialize_system_phase1(b, t);
00187     flat_initialize_system_phase2(b, b, eps, eta, dtout);
00188 
00189     // Initial output (to cerr, note) only if later output is anticipated.
00190 
00191     if (t_next <= t_end) {
00192         cerr << "Time = " << t << ",  steps = " << steps << "\n";
00193         print_energies(b, 0);
00194     }
00195 
00196     while (t <= t_end) {
00197         hdyn *bi;
00198         real ttmp = 1e100;
00199         bi = particle_to_move(b, ttmp);
00200 
00201         // Check for output and termination before taking the step.
00202 
00203         if (ttmp > t_next) {
00204             cerr << "Time = " << t << ",  steps = " << steps << "\n";
00205             print_energies(b, 1);
00206             t_next += dtout;
00207         }
00208 
00209         // Output a snapshot to cout at the scheduled time, or at end of run.
00210         // Use snap_cube_size to force all particles in the cube to have C.M.
00211         // at rest at the origin.
00212 
00213         if (ttmp > t_snap || ttmp > t_end) {
00214             if (snap_cube_size < VERY_LARGE_NUMBER)
00215                 shift_cm(b, snap_cube_size);
00216             put_node(cout, *b);
00217             cout << flush;
00218             t_snap += dt_snap;
00219         }
00220         if (ttmp > t_end)
00221             break;
00222 
00223         t = ttmp;
00224         predict_loworder_all(b, t);
00225         bi->clear_interaction();
00226         bi->flat_calculate_acc_and_jerk(b, eps * eps);
00227         bi->correct();
00228         bi->flat_update(eta, dtout);
00229         bi->store_old_force();
00230         steps++;
00231     }
00232 }
00233 
00234 void print_usage_and_exit()
00235 {
00236     cerr <<
00237          "usage: hermite -t # -a # -d # -D # -e # "
00238          << "[-c \"...\"]      "
00239          << "for t (time span),\n a (accuracy parameter ), "
00240          << "d (output interval), D (snapshot output interval), "
00241          << "e (softening length),\n";
00242     exit(1);
00243 }
00244 
00245 main(int argc, char **argv)
00246 {
00247     hdyn *b;                    // hdyn root node
00248 
00249     real delta_t = 10;          // time span of the integration
00250     real dtout = .25;           // output interval--make a power of 0.5
00251     real dt_snap;               // snap output interval
00252     real eps = 0.05;            // softening length               
00253     real eta = 0.05;            // time step parameter
00254 
00255     char *comment;              // comment string
00256 
00257     real snap_cube_size = VERY_LARGE_NUMBER;
00258 
00259     extern char *poptarg;
00260     int pgetopt(int, char **, char *), c;
00261 
00262     bool a_flag = FALSE;
00263     bool c_flag = FALSE;
00264     bool d_flag = FALSE;
00265     bool D_flag = FALSE;
00266     bool e_flag = FALSE;
00267     bool q_flag = FALSE;
00268     bool t_flag = FALSE;
00269 
00270     while ((c = pgetopt(argc, argv, "a:c:C:d:D:e:qt:")) != -1) {
00271         switch (c) {
00272             case 'a':
00273                 a_flag = TRUE;
00274                 eta = atof(poptarg);
00275                 break;
00276             case 'c':
00277                 c_flag = TRUE;
00278                 comment = poptarg;
00279                 break;
00280             case 'C':
00281                 snap_cube_size = atof(poptarg);
00282                 break;
00283             case 'd':
00284                 d_flag = TRUE;
00285                 dtout = atof(poptarg);
00286                 break;
00287             case 'D':
00288                 D_flag = TRUE;
00289                 dt_snap = atof(poptarg);
00290                 break;
00291             case 'e':
00292                 e_flag = TRUE;
00293                 eps = atof(poptarg);
00294                 break;
00295             case 'q':
00296                 q_flag = TRUE;
00297                 break;
00298             case 't':
00299                 t_flag = TRUE;
00300                 delta_t = atof(poptarg);
00301                 break;
00302             case '?':
00303                 print_usage_and_exit();
00304         }
00305     }
00306 
00307     if (!q_flag) {
00308 
00309         // Check input arguments and echo defaults.
00310 
00311         if (!t_flag)
00312             cerr << "default delta_t = " << delta_t << "\n";
00313         if (!d_flag)
00314             cerr << "default dtout = " << dtout << "\n";
00315         if (!a_flag)
00316             cerr << "default eta = " << eta << "\n";
00317         if (!e_flag)
00318             cerr << "default eps = " << eps << "\n";
00319     }
00320     if (!D_flag)
00321         dt_snap = delta_t;
00322 
00323     b = get_hdyn(cin);
00324 
00325     if (c_flag == TRUE)
00326         b->log_comment(comment);
00327     b->log_history(argc, argv);
00328 
00329     evolve_system(b, delta_t, eta, dtout, dt_snap, eps, snap_cube_size);
00330 }