nbody_sh1p.c

mpi编程

	   << dt_dia << ",\n  and snapshot output interval dt_out = "
	   << dt_out << "." << endl;

    real tcpu = MPI_Wtime();           // check CPU usage

    real epot = 0;                     // potential energy of the n-body system
    real coll_time = VERY_LARGE_NUMBER;// collision (close encounter) time scale
    get_acc_jerk_pot_coll(p, n, epot, coll_time, pipe);

    int nsteps = 0;               // number of integration time steps completed
    real einit;                   // initial total energy of the system

    write_diagnostics(p, n, t, epot, nsteps, einit,
		      true, x_flag, tcpu);
    if (init_out)                                    // flag for initial output
      put_snapshot(p, n, t, particletype);

    real t_dia = t + dt_dia;           // next time for diagnostics output
    real t_out = t + dt_out;           // next time for snapshot output
    real t_end = t + dt_tot;           // final time, to finish the integration

    while (true){
        while (t < t_dia && t < t_out && t < t_end){
            real dt_local = dt_param * coll_time;
	    real dt;
	    MPI_Allreduce(&dt_local, &dt, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); // synchronize time step
	    evolve_step(p, n, t, dt, epot, coll_time, pipe);
            nsteps++;
        }
        if (t >= t_dia){
            write_diagnostics(p, n, t, epot, nsteps,
                              einit, false, x_flag, tcpu);
            t_dia += dt_dia;
        }
        if (t >= t_out){
            put_snapshot(p, n, t, particletype);
            t_out += dt_out;
        }
        if (t >= t_end)
            break;
    }
}

/*-----------------------------------------------------------------------------
 *  evolve_step  --  takes one integration step for an N-body system, using the
 *                   Hermite algorithm.
 *-----------------------------------------------------------------------------
 */

void evolve_step(Particle p[], int n, real & t,
                 real dt, real & epot, real & coll_time,
		 void *pipe)
{

  Particle *po = new Particle[n];

    for (int i = 0; i < n ; i++)
        for (int k = 0; k < NDIM ; k++){
          po[i].pos[k] = p[i].pos[k];
          po[i].vel[k] = p[i].vel[k];
          po[i].acc[k] = p[i].acc[k];
          po[i].jerk[k] = p[i].jerk[k];
        }

    predict_step(p, n, dt);
    get_acc_jerk_pot_coll(p, n, epot, coll_time, pipe);
    correct_step(p, po, n, dt);
    t += dt;

    delete[] po;
}

/*-----------------------------------------------------------------------------
 *  predict_step  --  takes the first approximation of one Hermite integration
 *                    step, advancing the positions and velocities through a
 *                    Taylor series development up to the order of the jerks.
 *-----------------------------------------------------------------------------
 */

void predict_step(Particle p[], int n, real dt)
{
    for (int i = 0; i < n ; i++)
        for (int k = 0; k < NDIM ; k++){
            p[i].pos[k] += p[i].vel[k]*dt + p[i].acc[k]*dt*dt/2
                                      + p[i].jerk[k]*dt*dt*dt/6;
            p[i].vel[k] += p[i].acc[k]*dt + p[i].jerk[k]*dt*dt/2;
        }
}

/*-----------------------------------------------------------------------------
 *  correct_step  --  takes one iteration to improve the new values of position
 *                    and velocities, effectively by using a higher-order
 *                    Taylor series constructed from the terms up to jerk at
 *                    the beginning and the end of the time step.
 *-----------------------------------------------------------------------------
 */

void correct_step(Particle p[], Particle po[], int n, real dt)
{
    for (int i = 0; i < n ; i++)
        for (int k = 0; k < NDIM ; k++){
            p[i].vel[k] = po[i].vel[k] + (po[i].acc[k] + p[i].acc[k])*dt/2
                                       + (po[i].jerk[k] - p[i].jerk[k])*dt*dt/12;
            p[i].pos[k] = po[i].pos[k] + (po[i].vel[k] + p[i].vel[k])*dt/2
                                       + (po[i].acc[k] - p[i].acc[k])*dt*dt/12;
        }
}

/*-----------------------------------------------------------------------------
 *  get_acc_jerk_pot_coll  --  calculates accelerations and jerks, and as side
 *                             effects also calculates potential energy and
 *                             the time scale coll_time for significant changes
 *                             in local configurations to occur.
 *                                                  __                     __
 *                                                 |          -->  -->       |
 *               M                           M     |           r  . v        |
 *   -->          j    -->       -->          j    | -->        ji   ji -->  |
 *    a   ==  --------  r    ;    j   ==  -------- |  v   - 3 ---------  r   |
 *     ji     |-->  |3   ji        ji     |-->  |3 |   ji      |-->  |2   ji |
 *            | r   |                     | r   |  |           | r   |       |
 *            |  ji |                     |  ji |  |__         |  ji |     __|
 *                             
 *  note: it would be cleaner to calculate potential energy and collision time
 *        in a separate function.  However, the current function is by far the
 *        most time consuming part of the whole program, with a double loop
 *        over all particles that is executed every time step.  Splitting off
 *        some of the work to another function would significantly increase
 *        the total computer time (by an amount close to a factor two).
 *
 *  We determine the values of all four quantities of interest by walking
 *  through the system in a double {i,j} loop.  The first three, acceleration,
 *  jerk, and potential energy, are calculated by adding successive terms;
 *  the last, the estimate for the collision time, is found by determining the 
 *  minimum value over all particle pairs and over the two choices of collision
 *  time, position/velocity and sqrt(position/acceleration), where position and
 *  velocity indicate their relative values between the two particles, while
 *  acceleration indicates their pairwise acceleration.  At the start, the
 *  first three quantities are set to zero, to prepare for accumulation, while
 *  the last one is set to a very large number, to prepare for minimization.
 *       The integration loops only over half of the pairs, with j > i, since
 *  the contributions to the acceleration and jerk of particle j on particle i
 *  is the same as those of particle i on particle j, apart from a minus sign
 *  and a different mass factor.
 *-----------------------------------------------------------------------------
 */

void get_acc_jerk_pot_coll(Particle pl[], int nl, 
			   Particle po[], int no, 
			   real & epot, real & coll_time)
{

    real coll_time_q = VERY_LARGE_NUMBER;      // collision time to 4th power
    real coll_est_q;                           // collision time scale estimate
                                               // to 4th power (quartic)
    for (int i = 0; i < nl ; i++){
        for (int j = 0; j < no ; j++){            // rji[] is the vector from
	  if(pl[i].id!=po[j].id) {
            real rji[NDIM];                        // particle i to particle j
            real vji[NDIM];                        // vji[] = d rji[] / d t
            for (int k = 0; k < NDIM ; k++){
                rji[k] = po[j].pos[k] - pl[i].pos[k];
                vji[k] = po[j].vel[k] - pl[i].vel[k];
            }
            real r2 = 0;                           // | rji |^2
            real v2 = 0;                           // | vji |^2
            real rv_r2 = 0;                        // ( rij . vij ) / | rji |^2
            for (int k = 0; k < NDIM ; k++){
                r2 += rji[k] * rji[k];
                v2 += vji[k] * vji[k];
                rv_r2 += rji[k] * vji[k];
            }
            rv_r2 /= r2;
            real r = sqrt(r2);                     // | rji |
            real r3 = r * r2;                      // | rji |^3

// add the {i,j} contribution to the total potential energy for the system:

            epot -= pl[i].mass * po[j].mass / r;

// add the {j (i)} contribution to the {i (j)} values of acceleration and jerk:

            real da[3];                            // main terms in pairwise
            real dj[3];                            // acceleration and jerk
            for (int k = 0; k < NDIM ; k++){
                da[k] = rji[k] / r3;                           // see equations
                dj[k] = (vji[k] - 3 * rv_r2 * rji[k]) / r3;    // in the header
            }
            for (int k = 0; k < NDIM ; k++){
                pl[i].acc[k] += po[j].mass * da[k];           // using symmetry
		pl[i].jerk[k] += po[j].mass * dj[k];          // acceleration

		// in the original version pij = -pji for acc and jerk.
		// in this parallel version this is unpractical.
                //po[j].acc[k] -= pl[i].mass * da[k];      // find pairwise
		//po[j].jerk[k] -= pl[i].mass * dj[k];     // and jerk
            }

// first collision time estimate, based on unaccelerated linear motion:

            coll_est_q = (r2*r2) / (v2*v2);
            if (coll_time_q > coll_est_q)
                coll_time_q = coll_est_q;

// second collision time estimate, based on free fall:

            real da2 = 0;                                  // da2 becomes the 
            for (int k = 0; k < NDIM ; k++)                // square of the 
                da2 += da[k] * da[k];                      // pair-wise accel-
            double mij = pl[i].mass + po[j].mass;            // eration between
            da2 *= mij * mij;                              // particles i and j

            coll_est_q = r2/da2;
            if (coll_time_q > coll_est_q)
                coll_time_q = coll_est_q;
	  }
        }                                     
    }    
    // from q for quartic back to linear collision time and taking the minimum
    coll_time = min(coll_time, sqrt(sqrt(coll_time_q))); 
}                                             

/*-----------------------------------------------------------------------------
 *  get_acc_jerk_pot_coll  --  the distribution of particles over 
 *                             neighboring processors and subsequent force
 *                             calculation.
 *-----------------------------------------------------------------------------
 */

void get_acc_jerk_pot_coll(Particle p[], int n, real &epot, real &coll_time, 
			   void *pipe) {

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

    int rlen;
    Particle *recvbuf;

    for (int i = 0; i < n ; i++)
      for (int k = 0; k < NDIM ; k++) {
	p[i].acc[k] = p[i].jerk[k] = 0;
      }

    MPE_Pipe_start( pipe, p, n, 1 ); // load the initial sendbuffer 

    epot = 0;                        // initialize epot and coll_time
    coll_time = VERY_LARGE_NUMBER;
    get_acc_jerk_pot_coll(p, n, p, n, epot, coll_time);

    for (int step=1; step<size; step++) { // compute forces for other particles

      MPE_Pipe_push( pipe, (void**)&recvbuf, &rlen );        // get new data
      
      // Compute forces 
      get_acc_jerk_pot_coll(p, n, recvbuf, rlen, epot, coll_time); 
    }
}

/*-----------------------------------------------------------------------------
 *                                                                    \\   o
 *  end of file:  nbody_sh1.C                                         /\\'  O
 *                                                                   /\     |
 *=============================================================================
 */