flat_hermite.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00011 // flat_hermite.C: hermite scheme with simple (flat-tree) _dyn_ structure
00012 //
00013 // J. Makino    92-12-04:  Seems to work... Energy conservation looks good.
00014 // S. McMillan  93-04-02:  Cleaned up and modified command-line interface.
00015 // J. Makino    96-07-xx:  Seems not to work... Energy conservation looks bad.
00016 // S. McMillan  99-06-xx:  Seems to work... Energy conservation looks good.
00017 
00018 #include "_dyn_.h"
00019 
00020 local void accumulate_energies(_dyn_ * b,
00021                                real & epot, real & ekin, real & etot)
00022 {
00023     // Recursion is one level deep for flat trees...
00024 
00025     if (b->get_oldest_daughter()) {
00026 
00027         for_all_daughters(_dyn_, b, bb)
00028             accumulate_energies(bb, epot, ekin, etot);  // recursion...
00029 
00030     } else {
00031 
00032         real mi = b->get_mass();
00033         epot += 0.5 * mi * b->get_pot();
00034         vector vel = b->get_vel();
00035         ekin += 0.5 * mi * vel * vel;
00036 
00037     }
00038     etot = ekin + epot;
00039 }
00040 
00041 local void compute_energies(_dyn_ * b,
00042                             real & epot, real & ekin, real & etot)
00043 {
00044     epot = ekin = etot = 0;
00045     accumulate_energies(b, epot, ekin, etot);
00046 }
00047 
00048 local void print_energies(_dyn_ * b, int mode)
00049 {
00050     real epot, ekin, etot;
00051     static real old_etot = 0;
00052 
00053     compute_energies(b, epot, ekin, etot);
00054     cerr << "    Energies:  " << ekin << "  " << epot << "  " << etot;
00055 
00056     if (mode) {
00057         real de = (etot - old_etot) / etot;
00058         cerr << ";  de = " << de << endl;
00059     } else {
00060         old_etot = etot;
00061         cerr << endl;
00062     }
00063 }
00064 
00065 local _dyn_ *particle_to_move(_dyn_ * b, real & tmin)
00066 {
00067     // Scheduler returns a pointer to a single particle.  Since steps
00068     // are forced to be powers of 2 (see _dyn_ev.C), it would be better
00069     // if we returned a list of nodes to move, as in hdyn/evolve/kira.
00070 
00071     _dyn_ *bmin = NULL;
00072 
00073     if (b->get_oldest_daughter()) {
00074 
00075         for_all_daughters(_dyn_, b, bb) {
00076             real ttmp = VERY_LARGE_NUMBER;
00077             _dyn_ *btmp = particle_to_move(bb, ttmp);   // recursion...
00078             if (ttmp < tmin) {
00079                 tmin = ttmp;
00080                 bmin = btmp;
00081             }
00082         }
00083 
00084     } else {
00085 
00086         if (tmin > b->get_next_time()) {
00087             bmin = b;
00088             tmin = b->get_next_time();
00089         }
00090 
00091     }
00092 
00093     return bmin;
00094 }
00095 
00096 local void flat_initialize_system_phase1(_dyn_ * b,  real t)
00097 {
00098     if (b->get_oldest_daughter())
00099         for_all_daughters(_dyn_, b, d)
00100             flat_initialize_system_phase1(d, t);        // recursion...
00101 
00102     b->set_time(t);
00103     b->set_system_time(t);
00104     b->init_pred();
00105 }
00106 
00107 local void flat_initialize_system_phase2(_dyn_ * b, _dyn_ * root, real eps,
00108                                          real eta, real dtout)
00109 {
00110     if (b->get_oldest_daughter() != NULL) {
00111 
00112         for_all_daughters(_dyn_, b, bb)
00113             flat_initialize_system_phase2(bb, root,
00114                                           eps, eta, dtout); // recursion...
00115 
00116     } else {
00117 
00118         b->clear_interaction();
00119         b->flat_calculate_acc_and_jerk(root, eps * eps);
00120         b->store_old_force();
00121         b->flat_set_first_timestep(0.5 * eta, dtout);
00122 
00123     }
00124 }
00125 
00126 local void evolve_system(_dyn_ * b,     // root node
00127                          real delta_t,  // time span of the integration 
00128                          real eta,      // accuracy parameter 
00129                          real dtout,    // output time interval
00130                          real dt_snap,  // snapshot output interval
00131                          real eps,      // softening length             
00132                          real snap_cube_size)
00133 {
00134     real t = b->get_time();     // current time
00135     real t_end = t + delta_t;   // final time, at the end of the integration
00136     real t_next = t + dtout;    // time of next output
00137     real t_snap = t + dt_snap;  // time of next snapshot;
00138 
00139     int steps = 0;
00140 
00141     flat_initialize_system_phase1(b, t);
00142     flat_initialize_system_phase2(b, b, eps, eta, dtout);
00143 
00144     // Initial output (to cerr, note) only if later output is anticipated.
00145 
00146     if (t_next <= t_end) {
00147         cerr << "Time = " << t << ",  steps = " << steps << endl;
00148         print_energies(b, 0);
00149     }
00150 
00151     // Main integration loop.
00152 
00153     while (t <= t_end) {
00154 
00155         real ttmp = 1e100;
00156         _dyn_ *bi = particle_to_move(b, ttmp);
00157 
00158         // Check for output and termination before taking the step.
00159 
00160         if (ttmp > t_next) {
00161             cerr << "Time = " << t << ",  steps = " << steps << endl;
00162             print_energies(b, 1);
00163             t_next += dtout;
00164         }
00165 
00166         // Output a snapshot to cout at the scheduled time, or at end of run.
00167         // Use snap_cube_size to force all particles in the cube to have C.M.
00168         // at rest at the origin.
00169 
00170         if (ttmp > t_snap || ttmp > t_end) {
00171             put_node(cout, *b);
00172             cout << flush;
00173             t_snap += dt_snap;
00174         }
00175         if (ttmp > t_end) break;
00176 
00177         b->set_system_time(t=ttmp);
00178         predict_loworder_all(b, t);
00179 
00180         bi->clear_interaction();
00181 
00182         bi->flat_calculate_acc_and_jerk(b, eps * eps);
00183 
00184         bi->flat_correct();
00185         bi->flat_update(eta, dtout);
00186         bi->store_old_force();
00187 
00188         steps++;
00189     }
00190 }
00191 
00192 main(int argc, char **argv)
00193 {
00194     _dyn_ *b;                   // root node
00195 
00196     real delta_t = 10;          // time span of the integration
00197     real dtout = .25;           // output interval--make a power of 0.5
00198     real dt_snap;               // snap output interval
00199     real eps = 0.05;            // softening length               
00200     real eta = 0.05;            // time step parameter
00201 
00202     char *comment;              // comment string
00203 
00204     real snap_cube_size = VERY_LARGE_NUMBER;
00205 
00206     bool a_flag = FALSE;
00207     bool c_flag = FALSE;
00208     bool d_flag = FALSE;
00209     bool D_flag = FALSE;
00210     bool e_flag = FALSE;
00211     bool q_flag = FALSE;
00212     bool t_flag = FALSE;
00213 
00214     check_help();
00215 
00216     extern char *poptarg;
00217     int c;
00218     char* param_string = "a:c:C:d:D:e:qt:";
00219 
00220     while ((c = pgetopt(argc, argv, param_string)) != -1)
00221 
00222         switch (c) {
00223             case 'a':   a_flag = TRUE;
00224                         eta = atof(poptarg);
00225                         break;
00226             case 'c':   c_flag = TRUE;
00227                         comment = poptarg;
00228                         break;
00229             case 'C':   snap_cube_size = atof(poptarg);
00230                         break;
00231             case 'd':   d_flag = TRUE;
00232                         dtout = atof(poptarg);
00233                         break;
00234             case 'D':   D_flag = TRUE;
00235                         dt_snap = atof(poptarg);
00236                         break;
00237             case 'e':   e_flag = TRUE;
00238                         eps = atof(poptarg);
00239                         break;
00240             case 'q':   q_flag = TRUE;
00241                         break;
00242             case 't':   t_flag = TRUE;
00243                         delta_t = atof(poptarg);
00244                         break;
00245             case '?':   params_to_usage(cerr, argv[0], param_string);
00246                         get_help();
00247                         exit(1);
00248         }
00249 
00250     if (!q_flag) {
00251 
00252         // Check input arguments and echo defaults.
00253 
00254         if (!t_flag)
00255             cerr << "default delta_t = " << delta_t << endl;
00256         if (!d_flag)
00257             cerr << "default dtout = " << dtout << endl;
00258         if (!a_flag)
00259             cerr << "default eta = " << eta << endl;
00260         if (!e_flag)
00261             cerr << "default eps = " << eps << endl;
00262     }
00263     if (!D_flag)
00264         dt_snap = delta_t;
00265 
00266     b = get__dyn_(cin);
00267 
00268     if (c_flag == TRUE)
00269         b->log_comment(comment);
00270     b->log_history(argc, argv);
00271 
00272     evolve_system(b, delta_t, eta, dtout, dt_snap, eps, snap_cube_size);
00273 }