nbody_sh1p.c

mpi编程

               // 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 routines

using 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 << endl

const int NDIM = 3;                        // number of spatial dimensions
const real VERY_LARGE_NUMBER = 1e300;
const int root = 0;                        // identity of the root processor

// The Particle structure
typedef 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)]"