The following code evaluates accelerations on NI particles from other NJ particles by direct summation.
pt
int calculate_accel_by_grape6_separate_trial_noopen(int clusterid, int ni, double xi[][3], double vi[][3], int nj, double xj[][3], double vj[][3], double m[], double a[][3], double jerk[][3], double pot[], double eps2) { #define IPLIMIT MAXPIPELINESPERCHIP double ajtmp[3]; double jjtmp[3]; double j2jtmp[3]; double ti, tj, dtj; int j0, i0; int i,k,ii, iend; int flag[MAXPIPELINESPERCHIP+1]; int index[MAXPIPELINESPERCHIP+1]; double h2[MAXPIPELINESPERCHIP+1]; double eps2array[MAXPIPELINESPERCHIP+1]; int nharderror = 0; int ipmax = g6_npipes(); ti = 0.0;tj =0.0; dtj = 0.0078125; for(k=0;k<3;k++){ ajtmp[k] = 0.0; jjtmp[k] = 0.0; j2jtmp[k] = 0.0; } START: g6_set_ti(clusterid, 0.0); for(i=0;i<nj;i++){ g6_set_j_particle_mxonly(clusterid,i,i,m+i,xj[i]); } for(i=0;i<ipmax;i++){ h2[i] = eps2; eps2array[i] = eps2; } g6_flush_jp_buffer(clusterid); for(i=0;i<ni;i+=ipmax){ int error; iend = ipmax; if (iend+i > ni) iend = ni-i; for(ii=0;ii<iend;ii++)index[ii]=i+ii; g6calc_firsthalf(clusterid, nj, iend,index,&(xi[i]),&(vi[i]),&(a[i]),&(jerk[i]), pot+i,eps2, h2); if (error = g6calc_lasthalf(clusterid,nj,iend,index,&xi[i],&vi[i],eps2,h2, &a[i], &jerk[i], pot+i)){ nharderror ++; fprintf(stderr,"(calculate_accell_trial_noopen) hard error %d\n", error); g6_reinitialize(clusterid); if (nharderror < 10){ goto START; }else{ fprintf(stderr,"(calculate_accell_trial_noopen) too many errors... %d returning -1\n", nharderror); return -1; } } } return 0; }
The following is a C++ example to show some idea...
pt
g6_open(0); g6_initialize_jp_buffer(0,nbody); for(int i=0;i<nbody;i++)if (pb[i].get_grape_index() == -1) { pb[i].set_timestep(0.0078125); pb[i].set_grape_index(i); } vector j218 = vector(0.0,0.0,0.0); particle * bi = pb; for(int i = 0; i<nbody; i++){ vector j6 = bi->get_old_jerk()*ONE_SIXTH; vector a2 = bi->get_old_acc()*0.5L; g6_set_j_particle(grape6_id,bi->get_grape_index(), bi->get_index(), bi->get_time(), bi->get_timestep(), bi->get_grape_mass(), (real*)&j218, (real*)&j6, (real*)&a2, (real*)bi->pget_vel(), (real*)bi->pget_pos()); bi++; } }
void calculate_acc_and_jerk_for_list_on_grape6(particle_system* ps, particle* pb, int nbody, int nbh, node_time* nt, int n_next, real eps2) { int error; do{ error = 0; if (grape6_open == 0) { // actual initialization should be performed here... } real sys_t = nt[0].next_time; for (int i = 0; i < n_next; i++) { particle *bi = nt[i].pptr; bi->predict_loworder(sys_t); } g6_set_ti(grape6_id,sys_t); for (int i = 0; i < n_next; i+= npipes) { particle *bi = nt[i].pptr; int np = npipes; if (i+np > n_next) np = n_next - i; for (int ii = 0; ii < np; ii++){ particle *bi = nt[i+ii].pptr; gindex[ii] = bi->get_index(); xi[ii] = bi->get_pred_pos(); veli[ii] = bi->get_pred_vel(); acci[ii] = bi->get_acc(); jerki[ii] = bi->get_jerk(); poti[ii] = bi->get_pot(); } if ( (np > 0) && (error == 0)) g6calc_firsthalf(grape6_id,nbody, np, gindex, cvector xi, cvector veli, cvector acci, cvector jerki, poti, eps2, h2i); if ( (np > 0) && (error == 0)){ error = g6calc_lasthalf(grape6_id, nbody, np, gindex, cvector xi, cvector veli, eps2, h2i, cvector acci, cvector jerki, poti); if (error){ error = 1; }else{ for (int ii = 0; ii < np; ii++){ particle *bi = nt[i+ii].pptr; bi->set_acc(acci[ii]); bi->set_jerk(jerki[ii]); bi->set_pot(poti[ii]); } } } } if (error){ cerr << "GRAPE hardware error" <<endl; grape6_reinitialize(grape6_id); grape6_open = 0; } }while (error); }
pt
void particle_system::update_grape_data(int n_next) { vector j218 = vector(0.0,0.0,0.0); for (int i = 0; i < n_next; i++) { particle *bi = nt[i].pptr; vector j6 = bi->get_jerk()*ONE_SIXTH; vector a2 = bi->get_acc()*0.5L; g6_set_j_particle(grape6_id,bi->get_grape_index(), bi->get_index(), bi->get_time(), bi->get_timestep(), bi->get_grape_mass(), (real*)&j218, (real*)&j6, (real*)&a2, (real*)bi->pget_vel(), (real*)bi->pget_pos()); } g6_flush_jp_buffer(grape6_id); }