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);
}