nbody_sh1p.c

mpi编程

                           << " [-d step_size_control_parameter]\n"
                           << "         [-e diagnostics_interval]"
                           << " [-o output_interval]\n"
                           << "         [-t total_duration]"
                           << " [-i (start output at t = 0)]\n"
                           << "         [-x (extra debugging diagnostics)]"
                           << endl;
                      return false;         // execution should stop after help
            case 'd': dt_param = atof(optarg);
                      break;
            case 'e': dt_dia = atof(optarg);
                      break;
            case 'i': i_flag = true;
                      break;
            case 'o': dt_out = atof(optarg);
                      break;
            case 't': dt_tot = atof(optarg);
                      break;
            case 'x': x_flag = true;
                      break;
            case '?': cerr << "usage: " << argv[0]
                           << " [-h (for help)]"
                           << " [-d step_size_control_parameter]\n"
                           << "         [-e diagnostics_interval]"
                           << " [-o output_interval]\n"
                           << "         [-t total_duration]"
                           << " [-i (start output at t = 0)]\n"
                           << "         [-x (extra debugging diagnostics)]"
                           << endl;
                      return false;        // execution should stop after error
            }

    return true;                         // ready to continue program execution
}

/*-----------------------------------------------------------------------------
 *  get_snapshot  --  reads a single snapshot with the root processor 
 *                    from the input stream cin.
 *
 *  note: in this implementation, particle number, time and the data for
 *        individual particles are read in.
 *        All communication goes via the root (0) process.
 *  note: The routine also declared the MPI_Datatype particletype which is
 *        the contianer for the particle data.
 *        This structure is used for further MPI comminication.
 *-----------------------------------------------------------------------------
 */

Particle *get_snapshot(int &n, real &t,
		       MPI_Datatype &particletype)
{

  int rank, size;
  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
  MPI_Comm_size( MPI_COMM_WORLD, &size );

  Particle *p_tmp;
  if(rank==root) {
    cerr << "Reading snapshot" << endl;
    cin >> n;
    cin >> t;

    p_tmp = new Particle[n];

    for(int i=0; i<n; i++) {
      p_tmp[i].id = i;
      cin >> p_tmp[i].mass;                         // mass of particle i
      for (int k = 0; k < NDIM; k++)
	cin >> p_tmp[i].pos[k];                 // position of particle i
      for (int k = 0; k < NDIM; k++)
	cin >> p_tmp[i].vel[k];                 // velocity of particle i
    }
  }

  MPI_Bcast(&n,1,MPI_INT,root,MPI_COMM_WORLD); // broadcasts particle number

  int n_local = (int)(floor(1.0*n/size));
  if(n != n_local*size && rank==root) {
    cerr << "WARNING: Paticle number in input is not a mulitple of the number of processors." 
	 << endl;
    cerr << "         Action: Reduce particle number to n = " 
	 << n_local*size << "." << endl;
  }
  Particle *p = new Particle[n_local];

  // defining the particletype
  int inputblockcounts[2] = {1, 13};
  MPI_Datatype ntypes[] = {MPI::INT, MPI::DOUBLE};
  MPI::Aint displs[2];
  MPI_Address(&p[0].id, &displs[0]);
  MPI_Address(&p[0].mass, &displs[1]);
  displs[1] -= displs[0];        //make them relative
  displs[0] = 0;

  MPI_Type_struct(2, inputblockcounts, displs, ntypes, &particletype);
  MPI_Type_commit(&particletype);

  // Distribute the particles over the processors
  MPI_Scatter(p_tmp,n_local,particletype,p,n_local,particletype,root,MPI_COMM_WORLD);
  n = n_local;
  MPI_Bcast(&t,1,MPI_DOUBLE,root,MPI_COMM_WORLD);

  delete []p_tmp;
  return p;
}

/*-----------------------------------------------------------------------------
 *  put_snapshot  --  writes a single snapshot on the output stream cout.
 *
 *  note: Communication goes via the root (0) proceesor
 *-----------------------------------------------------------------------------
 */

void put_snapshot(Particle p[], int n, real t,
		  MPI_Datatype particletype)
{

    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    int ntot;
    MPI_Allreduce(&n, &ntot, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );

    Particle *p_all;
    if(rank==root) {
      p_all = new Particle[ntot];
    }
    MPI_Gather(p,n,particletype,p_all,n,particletype,root,MPI_COMM_WORLD);

    cout.precision(16);                       // full double precision
    if(rank==root) {
      cout << ntot << endl;                     // N, total particle number
      cout << t << endl;                        // current time
    for (int i = 0; i < ntot ; i++){
      cout << p_all[i].mass;                      // mass of particle i
        for (int k = 0; k < NDIM; k++)
            cout << ' ' << p_all[i].pos[k];         // position of particle i
        for (int k = 0; k < NDIM; k++)
            cout << ' ' << p_all[i].vel[k];         // velocity of particle i
        cout << endl;
    }
    delete []p_all;
    }
}

    
/*-----------------------------------------------------------------------------
 *  write_diagnostics  --  writes diagnostics on the error stream cerr:
 *                         current time; number of integration steps so far;
 *                         kinetic, potential, and total energy; absolute and
 *                         relative energy errors since the start of the run.
 *                         If x_flag (x for eXtra data) is true, all internal
 *                         data are dumped for each particle (mass, position,
 *                         velocity, acceleration, and jerk).
 *
 *  note: the kinetic energy is calculated here, while the potential energy is
 *        calculated in the function get_acc_jerk_pot_coll() and broadcasted
 *        by the other processors.
 *
 *  note: Communication goes via the root (0) proceesor
 *-----------------------------------------------------------------------------
 */

void write_diagnostics(Particle p[], int n, real t, real epot_local,
                       int nsteps, real & einit, bool init_flag,
                       bool x_flag, real &tcpu)
{

    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    real ekin_local = 0;       // kinetic energy of the n-body system
    for (int i = 0; i < n ; i++)
        for (int k = 0; k < NDIM ; k++)
            ekin_local += 0.5 * p[i].mass * p[i].vel[k] * p[i].vel[k];

    real ekin;
    MPI_Allreduce(&ekin_local, &ekin, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );

    real epot;
    MPI_Allreduce(&epot_local, &epot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
    epot *= 0.5;  // against double counting

    real etot = ekin + epot;             // total energy of the n-body system

    if (init_flag)                       // at first pass, pass the initial
        einit = etot;                    // energy back to the calling function

    tcpu = MPI_Wtime() - tcpu;
    if(rank==0) {


    cerr << "at time t = " << t << ", after " << nsteps
         << " steps (CPU = " << tcpu << "): \n  E_kin = " << ekin
         << " , E_pot = " << epot
         << " , E_tot = " << etot << endl;
    cerr << "                "
         << "absolute energy error: E_tot - E_init = "
         << etot - einit << endl;
    cerr << "                "
         << "relative energy error: (E_tot - E_init) / E_init = "
         << (etot - einit) / einit << endl;
    }
    if (x_flag){
        cerr << "  for debugging purposes, here is the internal data "
             << "representation:\n";
        for (int i = 0; i < n ; i++){
            cerr << "    internal data for particle " << i+1 << " : " << endl;
            cerr << "      ";
            cerr << p[i].id << " " << p[i].mass;
            for (int k = 0; k < NDIM; k++)
                cerr << ' ' << p[i].pos[k];
            for (int k = 0; k < NDIM; k++)
                cerr << ' ' << p[i].vel[k];
            for (int k = 0; k < NDIM; k++)
                cerr << ' ' << p[i].acc[k];
            for (int k = 0; k < NDIM; k++)
                cerr << ' ' << p[i].jerk[k];
            cerr << endl;
        }
    }
}
 
/*-----------------------------------------------------------------------------
 *  evolve  --  integrates an N-body system, for a total duration dt_tot.
 *              Snapshots are sent to the standard output stream once every
 *              time interval dt_out.  Diagnostics are sent to the standard
 *              error stream once every time interval dt_dia.
 *
 *  note: the integration time step, shared by all particles at any given time,
 *        is variable.  Before each integration step we use coll_time (short
 *        for collision time, an estimate of the time scale for any significant
 *        change in configuration to happen), multiplying it by dt_param (the
 *        accuracy parameter governing the size of dt in units of coll_time),
 *        to obtain the new time step size.
 *
 *  Before moving any particles, we start with an initial diagnostics output
 *  and snapshot output if desired.  In order to write the diagnostics, we
 *  first have to calculate the potential energy, with get_acc_jerk_pot_coll().
 *  That function also calculates accelerations, jerks, and an estimate for the
 *  collision time scale, all of which are needed before we can enter the main
 *  integration loop below.
 *       In the main loop, we take as many integration time steps as needed to
 *  reach the next output time, do the output required, and continue taking
 *  integration steps and invoking output this way until the final time is
 *  reached, which triggers a `break' to jump out of the infinite loop set up
 *  with `while(true)'.
 *-----------------------------------------------------------------------------
 */

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)
{
    int rank, size;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    MPI_Comm_size( MPI_COMM_WORLD, &size );

    if(rank==root) 
      cerr << "Starting a Hermite integration for a " << n*size
	   << "-body system,\n  from time t = " << t 
	   << " with time step control parameter dt_param = " << dt_param
	   << "  until time " << t + dt_tot 
	   << " ,\n  with diagnostics output interval dt_dia = "