next up previous contents
Next: 7 Error recovery Up: grape6user Previous: 5 Reference Manual for

Subsections


6 Example

6.1 Shared timestep

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

6.2 Block timestep

The following is a C++ example to show some idea...

6.2.1 Initialization

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

6.2.2 Force calculation

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

6.3 Update memory

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



Jun Makino
平成17年1月31日