// Time-stamp: <2002-01-18 21:51:36 piet>              // Time-stamp: <2003-21-10 03:22:41 spz>             //================================================================            //                                                                |           //           /__----__                         ........            |          //       .           \                     ....:        :.          |         //       :                 _\|/_         ..:                         |        //       :                   /|\         :                     _\|/_  |       //  ___   ___                  _____                      ___    /|\   |      //  /     |   \    /\ \     / |   |   \  / |        /\    |   \         |     //  |   __ |___/   /  \ \   /  |   |    \/  |       /  \   |___/         |    //   |    | |  \   /____\ \ /   |   |    /   |      /____\  |   \     \/  |   //     \___| |   \ /      \ V    |   |   /    |____ /      \ |___/     |   |  //                                                                      /   | //              :                       _/|     :..                    |/    |//                :..               ____/           :....          ..         |/*   o   //          :.    _\|/_     /                   :........:           | *  O  `//\                 /|\                                               | *  |     /\                                                                  | *============================================================================= * *  nbody_sh1p.C: an N-body integrator with a shared but variable time step *                (the same for all particles but changing in time), using *                the Hermite integration scheme. *                         *                ref.: Hut, P., Makino, J. & McMillan, S., 1995, *                      Astrophysical Journal Letters 443, L93-L96. *                 *                The original version (nbody_sh1.C) has been adapted  *                for a parallel ring algorithm using the MPI library.  *_____________________________________________________________________________ * *  usage: nbody_sh1p [-h (for help)] [-d step_size_control_parameter] *                    [-e diagnostics_interval] [-o output_interval] *                    [-t total_duration] [-i (start output at t = 0)] *                    [-x (extra debugging diagnostics)] *  *         "step_size_control_parameter" is a coefficient determining the *            the size of the shared but variable time step for all particles * *         "diagnostics_interval" is the time between output of diagnostics, *            in the form of kinetic, potential, and total energy; with the *            -x option, a dump of the internal particle data is made as well *  *         "output_interval" is the time between successive snapshot outputs * *         "total_duration" is the integration time, until the program stops * *         Input and output are written from the standard i/o streams.  Since *         all options have sensible defaults, the simplest way to run the code *         is by only specifying the i/o files for the N-body snapshots: * *            nbody_sh1 < data.in > data.out * *         The diagnostics information will then appear on the screen. *         To capture the diagnostics information in a file, capture the *         standard error stream as follows: * *            (nbody_sh1 < data.in > data.out) >& data.err * *  Note: if any of the times specified in the -e, -o, or -t options are not an *        an integer multiple of "step", output will occur slightly later than *        predicted, after a full time step has been taken.  And even if they *        are integer multiples, round-off error may induce one extra step. *_____________________________________________________________________________ * *  External data format: * *     The program expects input of a single snapshot of an N-body system, *     in the following format: the number of particles in the snapshot n; *     the time t; mass mi, position ri and velocity vi for each particle i, *     with position and velocity given through their three Cartesian *     coordinates, divided over separate lines as follows: * *                      n *                      t *                      m1 r1_x r1_y r1_z v1_x v1_y v1_z *                      m2 r2_x r2_y r2_z v2_x v2_y v2_z *                      ... *                      mn rn_x rn_y rn_z vn_x vn_y vn_z * *     Output of each snapshot is written according to the same format. * *  Internal data format: * *     The data for an N-body system is stored internally as a 1-dimensional *     array for particle characteristics collected in a structure. *     The structure contains particle identity, mass, and 1-dimensional  *     arrays for the position, velocity, acceleration and jerk. *     The particle identity is added in the MPI implemetation to prevent  *     computing acceleration and force on itself. *_____________________________________________________________________________ * *  Extera about the MPI implementation * *     The topology assumes a 1-dimensional ring of processors with  *     int(N/p) particles per processor. At the moment only p*int(N/p) *     particles can be integrated; excess of p*int(N/p) particles are ignored  *     in the input and particle residtribution. *     All I/O is via the root processor (0).  * *     Possible adjustments and improvement: *     - It is assumed that each star has a unique identifier, which is  *       provided at runtime. *     - For the comminucation all particle information is passed to other  *       processes, where this is not always required. *     - Since each processor works with N/p particles, little memory per  *       node is requires, the root node, however, requires to store an  *       entire snapshot at I/O. *     - Acceleration and jerk are computed for all particles individually,  *       the symmetry in the force calculation is not used. *_____________________________________________________________________________ * *  Compiling * *    Successfull compilation requires the presence of MPI or some compatible  *    derivative such as MPICH. *    Here we give the command line example for compilation with MPICH  *    assuming that it is installed on your system on /usr/local/MPI * *    First: create the pipe.o object file: *    %> /usr/local/MPI/bin/mpiCC -g -O2 -I/usr/local/MPI/include -c pipe.C * *    Second: compile nbody_sh1p.C and link it with pipe.o: *    %> /usr/local/MPI/bin/mpiCC -g -O2 -I/usr/local/MPI/include \ *          nbody_sh1p.C -L/usr/local/MPI/lib -lmpich pipe.o -o nbody_sh1p *_____________________________________________________________________________ * *  Testrun *  *    Test run with 6 processors by the following command line *    %> /usr/local/MPI/bin/mpirun -np 6 ./nbody_sh1p < n24body.in *_____________________________________________________________________________ * *    version 1:  Jan 2002   Piet Hut, Jun Makino *    version 2:  Okt 2003   Simon Portegies Zwart *_____________________________________________________________________________ */#include  <iostream>#include  <cmath>                          // to include sqrt(), etc.#include  <cstdlib>                        // for atoi() and atof()#include  <unistd.h>                       // for getopt()#include  <mpi.h>                          // include the MPI library#include  "pipe.h"                         // the MPI pipe routinesusing namespace std;typedef double  real;                      // "real" as a general name for the                                           // standard floating-point data type// practical for debugging#define PRI(x) {for (int __pri__ = 0; __pri__ < x; __pri__++) cerr << " ";}#define PR(x)  cerr << #x << " = " << x << " "#define PRC(x) cerr << #x << " = " << x << ",  "#define PRL(x) cerr << #x << " = " << x << endlconst int NDIM = 3;                        // number of spatial dimensionsconst real VERY_LARGE_NUMBER = 1e300;const int root = 0;                        // identity of the root processor// The Particle structuretypedef struct {  int id;  real mass;  real pos[NDIM];  real vel[NDIM];  real acc[NDIM];  real jerk[NDIM];} Particle;void correct_step(Particle p[], Particle po[], int n, real dt);void evolve(Particle p[],            int n, real & t, real dt_param, real dt_dia, real dt_out,            real dt_tot, bool init_out, bool x_flag, void *pipe,	    MPI_Datatype particletype);void evolve_step(Particle p[], int n, real & t,                 real dt, real & epot, real & coll_time,		 void *pipe);void get_acc_jerk_pot_coll(Particle pl[], int nl, 			   Particle po[], int no, 			   real & epot, real & coll_time);Particle * get_snapshot(int &n, real &t, MPI_Datatype &particletype);void predict_step(Particle p[], int n, real dt);void put_snapshot(Particle p[], int n, real t, MPI_Datatype particletype);bool read_options(int argc, char *argv[], real & dt_param, real & dt_dia,                  real & dt_out, real & dt_tot, bool & i_flag, bool & x_flag);void write_diagnostics(Particle p[], int n, real t, real epot,                       int nsteps, real & einit, bool init_flag,                       bool x_flag, real &tcpu);#define loop(idx,last) for (idx = 0; idx < last ; idx++)void uniform(real a, real *x);void uniform(Particle p[], int n);void get_acc_jerk_pot_coll(Particle p[], int n, real &epot, real &coll_time, 			   void *pipe);/*----------------------------------------------------------------------------- *  main  --  reads option values, reads a snapshot, and launches the *            integrator *----------------------------------------------------------------------------- */int main(int argc, char *argv[]){    // initialize MPI    int rank, size;    MPI_Init( &argc, &argv );    MPI_Comm_rank( MPI_COMM_WORLD, &rank );    MPI_Comm_size( MPI_COMM_WORLD, &size );    real  dt_param = 0.03;     // control parameter to determine time step size    real  dt_dia = 1;          // time interval between diagnostics output    real  dt_out = 1;          // time interval between output of snapshots    real  dt_tot = 10;         // duration of the integration    bool  init_out = false;    // if true: snapshot output with start at t = 0                               //          with an echo of the input snapshot    bool  x_flag = false;      // if true: extra debugging diagnostics output    if (! read_options(argc, argv, dt_param, dt_dia, dt_out, dt_tot, init_out,                       x_flag))        return 1;                // halt criterion detected by read_options()    int n;    real t;    MPI_Datatype particletype;    Particle *p = get_snapshot(n, t, particletype);    real noutp = 1;    real dt;    put_snapshot(p, n, t, particletype);    void * pipe;  // Create a MPI pipe for a 1-dimensional ring topology    MPE_Pipe_create( MPI_COMM_WORLD, particletype, n, &pipe );    evolve(p, n, t, dt_param, dt_dia, dt_out, dt_tot, init_out,	   x_flag, pipe, particletype);    delete []p;    MPE_Pipe_free( &pipe );             // Clean up MPI    MPI_Type_free( &particletype );    MPI_Finalize();}/*----------------------------------------------------------------------------- *  read_options  --  reads the command line options, and implements them. * *  note: when the help option -h is invoked, the return value is set to false, *        to prevent further execution of the main program; similarly, if an *        unknown option is used, the return value is set to false. *----------------------------------------------------------------------------- */bool read_options(int argc, char *argv[], real & dt_param, real & dt_dia,                  real & dt_out, real & dt_tot, bool & i_flag, bool & x_flag){    int c;    while ((c = getopt(argc, argv, "hd:e:o:t:ix")) != -1)        switch(c){            case 'h': cerr << "usage: " << argv[0]                           << " [-h (for help)]"