Robust/src/Benchmarks/Prefetch/Moldyn/javasingle/JGFMolDynBench.java

   1 /**************************************************************************
   2  *                                                                         *
   3  *         Java Grande Forum Benchmark Suite - Thread Version 1.0          *
   4  *                                                                         *
   5  *                            produced by                                  *
   6  *                                                                         *
   7  *                  Java Grande Benchmarking Project                       *
   8  *                                                                         *
   9  *                                at                                       *
  10  *                                                                         *
  11  *                Edinburgh Parallel Computing Centre                      *
  12  *                                                                         *
  13  *                email: epcc-javagrande@epcc.ed.ac.uk                     *
  14  *                                                                         *
  15  *                                                                         *
  16  *      This version copyright (c) The University of Edinburgh, 2001.      *
  17  *                         All rights reserved.                            *
  18  *                                                                         *
  19  **************************************************************************/
  20 public class JGFMolDynBench {
  21   public int ITERS;
  22   public double LENGTH;
  23   public double m;
  24   public double mu;
  25   public double kb;
  26   public double TSIM;
  27   public double deltat;
  28
  29   public int PARTSIZE;
  30
  31   public DoubleWrapper[] epot;
  32   public DoubleWrapper[] vir;
  33   public DoubleWrapper[] ek;
  34
  35   int size,mm;
  36   int[] datasizes;
  37
  38   public int interactions;
  39   public IntWrapper[] interacts;
  40
  41   public int nthreads;
  42   public JGFInstrumentor instr;
  43
  44   public JGFMolDynBench(int nthreads) {
  45     this.nthreads=nthreads;
  46   }
  47
  48   public void JGFsetsize(int size){
  49     this.size = size;
  50   }
  51
  52   public void JGFinitialise(){
  53     interactions = 0;
  54     datasizes = new int[3];
  55     datasizes[0] = 8;
  56     datasizes[1] = 13;
  57     datasizes[2] = 18;
  58
  59     mm = datasizes[size];
  60     PARTSIZE = mm*mm*mm*4;
  61     ITERS = 100;
  62     LENGTH = 50e-10;
  63     m = 4.0026;
  64     mu = 1.66056e-27;
  65     kb = 1.38066e-23;
  66     TSIM = 50;
  67     deltat = 5e-16;
  68   }
  69
  70   public static void JGFapplication(JGFMolDynBench mold) {
  71     // Create new arrays
  72     //BarrierServer mybarr;
  73     /*
  74        int[] mid = new int[8];
  75        mid[0] = (128<<24)|(195<<16)|(175<<8)|84;//dw-10
  76        mid[1] = (128<<24)|(195<<16)|(175<<8)|85;//dw-11
  77        mid[2] = (128<<24)|(195<<16)|(175<<8)|86;//dw-12
  78        mid[3] = (128<<24)|(195<<16)|(175<<8)|87;//dw-13
  79        mid[4] = (128<<24)|(195<<16)|(175<<8)|88;//dw-14
  80        mid[5] = (128<<24)|(195<<16)|(175<<8)|89;//dw-15
  81        mid[6] = (128<<24)|(195<<16)|(175<<8)|90;//dw-16
  82        mid[7] = (128<<24)|(195<<16)|(175<<8)|91;//dw-17
  83        */
  84
  85     double sh_force [][];
  86     double sh_force2 [][][];
  87     int partsize, numthreads;
  88     partsize = mold.PARTSIZE;
  89     numthreads = mold.nthreads;
  90     //mybarr = new BarrierServer(numthreads);
  91     //mybarr.start(mid[0]);
  92
  93     sh_force = new double[3][partsize];
  94     sh_force2 = new double[3][numthreads][partsize];
  95     mold.epot = new DoubleWrapper[numthreads];
  96     mold.vir  = new DoubleWrapper[numthreads];
  97     mold.ek   = new DoubleWrapper[numthreads];
  98     mold.interacts = new IntWrapper[numthreads];
  99     for(int i=0;i<numthreads;i++) {
 100       mold.epot[i]=new DoubleWrapper();
 101       mold.vir[i]=new DoubleWrapper();
 102       mold.ek[i]=new DoubleWrapper();
 103       mold.interacts[i]=new IntWrapper();
 104     }
 105
 106     // spawn threads
 107     MDWrap[] thobjects = new MDWrap[numthreads];
 108
 109     for(int i=0;i<numthreads;i++) {
 110       thobjects[i] = new MDWrap(new mdRunner(i,mold.mm,sh_force,sh_force2,mold.nthreads,mold));
 111     }
 112
 113     /*
 114        boolean waitfordone=true;
 115        while(waitfordone) {
 116        if (mybarr.done)
 117        waitfordone=false;
 118        }
 119        */
 120
 121     for(int i=0;i<numthreads;i++) {
 122       //thobjects[i].md.start(mid[i]);
 123       thobjects[i].md.run();
 124     }
 125   }
 126
 127   public void JGFvalidate(){
 128     double[] refval = new double[2];
 129     refval[0] = 1731.4306625334357;
 130     refval[1] = 7397.392307839352;
 131     double dev = Math.fabs(ek[0].d - refval[size]);
 132     if (dev > 1.0e-10 ){
 133       //System.printString("Validation failed\n");
 134       //System.printString("Kinetic Energy = " + (long)ek[0] + "  " + (long)dev + "  " + size + "\n");
 135     }
 136   }
 137 }
 138
 139 class mdRunner {
 140
 141   double count;
 142   int id,i,j,k,lg,mm;
 143   double l,rcoff,rcoffs,side,sideh,hsq,hsq2,vel,velt;
 144   double a,r,sum,tscale,sc,ekin,ts,sp;
 145   double den;
 146   double tref;
 147   double h;
 148   double vaver,vaverh,rand;
 149   double etot,temp,pres,rp;
 150   double u1, u2, s, xx, yy, zz;
 151   double xvelocity, yvelocity, zvelocity;
 152
 153   double [][] sh_force;
 154   double [][][] sh_force2;
 155
 156   int ijk,npartm,iseed,tint;
 157   int irep;
 158   int istop;
 159   int iprint;
 160
 161   JGFMolDynBench mymd;
 162   int nthreads;
 163
 164   public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,
 165       int nthreads, JGFMolDynBench mymd) {
 166     this.id=id;
 167     this.mm=mm;
 168     this.sh_force=sh_force;
 169     this.sh_force2=sh_force2;
 170     this.nthreads = nthreads;
 171     this.mymd = mymd;
 172     count = 0.0;
 173     den = 0.83134;
 174     tref = 0.722;
 175     h = 0.064;
 176     irep = 10;
 177     istop = 19;
 178     iprint = 10;
 179   }
 180
 181   public void init(particle[] one, int mdsize) {
 182     int id=this.id;
 183     for (int lg=0; lg<=1; lg++) {
 184       for (int i=0; i<mm; i++) {
 185         for (int j=0; j<mm; j++) {
 186           for (int k=0; k<mm; k++) {
 187             one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
 188                 xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
 189             ijk = ijk + 1;
 190           }
 191         }
 192       }
 193     }
 194
 195     for (int lg=1; lg<=2; lg++) {
 196       for (int i=0; i<mm; i++) {
 197         for (int j=0; j<mm; j++) {
 198           for (int k=0; k<mm; k++) {
 199             one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
 200                 (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
 201             ijk = ijk + 1;
 202           }
 203         }
 204       }
 205     }
 206
 207     /* Initialise velocities */
 208
 209     iseed = 0;
 210     double v1 = 0.0;
 211     double v2 = 0.0;
 212     random randnum = new random(iseed,v1,v2);
 213
 214
 215     for (int i=0; i<mdsize; i+=2) {
 216       r  = randnum.seed();
 217       one[i].xvelocity = r*randnum.v1;
 218       one[i+1].xvelocity  = r*randnum.v2;
 219     }
 220
 221     for (int i=0; i<mdsize; i+=2) {
 222       r  = randnum.seed();
 223       one[i].yvelocity = r*randnum.v1;
 224       one[i+1].yvelocity  = r*randnum.v2;
 225     }
 226
 227     for (int i=0; i<mdsize; i+=2) {
 228       r  = randnum.seed();
 229       one[i].zvelocity = r*randnum.v1;
 230       one[i+1].zvelocity  = r*randnum.v2;
 231     }
 232
 233
 234     /* velocity scaling */
 235
 236     ekin = 0.0;
 237     sp = 0.0;
 238
 239     for(int i=0;i<mdsize;i++) {
 240       sp = sp + one[i].xvelocity;
 241     }
 242     sp = sp / mdsize;
 243
 244     for(int i=0;i<mdsize;i++) {
 245       one[i].xvelocity = one[i].xvelocity - sp;
 246       ekin = ekin + one[i].xvelocity*one[i].xvelocity;
 247     }
 248
 249     sp = 0.0;
 250     for(int i=0;i<mdsize;i++) {
 251       sp = sp + one[i].yvelocity;
 252     }
 253     sp = sp / mdsize;
 254
 255     for(int i=0;i<mdsize;i++) {
 256       one[i].yvelocity = one[i].yvelocity - sp;
 257       ekin = ekin + one[i].yvelocity*one[i].yvelocity;
 258     }
 259
 260
 261     sp = 0.0;
 262     for(int i=0;i<mdsize;i++) {
 263       sp = sp + one[i].zvelocity;
 264     }
 265     sp = sp / mdsize;
 266
 267     for(int i=0;i<mdsize;i++) {
 268       one[i].zvelocity = one[i].zvelocity - sp;
 269       ekin = ekin + one[i].zvelocity*one[i].zvelocity;
 270     }
 271
 272     ts = tscale * ekin;
 273     sc = h * Math.sqrt(tref/ts);
 274
 275
 276     for(int i=0;i<mdsize;i++) {
 277
 278       one[i].xvelocity = one[i].xvelocity * sc;
 279       one[i].yvelocity = one[i].yvelocity * sc;
 280       one[i].zvelocity = one[i].zvelocity * sc;
 281
 282     }
 283
 284   }
 285
 286   public void doinit(int mdsize) {
 287     for(int j=0;j<3;j++) {
 288       double[] sh=sh_force[j];
 289       for (int i=0;i<mdsize;i++) {
 290         sh[i] = 0.0;
 291       }
 292     }
 293   }
 294
 295
 296   public void doinit2(int mdsize) {
 297     for(int k=0;k<3;k++) {
 298       double[] sh=sh_force[k];
 299       double [][] sha=sh_force2[k];
 300       for(int j=0;j<nthreads;j++) {
 301         double[] sha2=sha[j];
 302         for(int i=0;i<mdsize;i++) {
 303           sh[i] += sha2[i];
 304         }
 305       }
 306     }
 307
 308     for(int k=0;k<3;k++) {
 309       double [][] sh1=sh_force2[k];
 310       for(int j=0;j<nthreads;j++) {
 311         double[] sh2=sh1[j];
 312         for(int i=0;i<mdsize;i++) {
 313
 314
 315           sh2[i] = 0.0;
 316         }
 317       }
 318     }
 319
 320     for(int j=1;j<nthreads;j++) {
 321       mymd.epot[0].d += mymd.epot[j].d;
 322       mymd.vir[0].d += mymd.vir[j].d;
 323     }
 324     for(int j=1;j<nthreads;j++) {
 325       mymd.epot[j].d = mymd.epot[0].d;
 326       mymd.vir[j].d = mymd.vir[0].d;
 327     }
 328     for(int j=0;j<nthreads;j++) {
 329       mymd.interactions += mymd.interacts[j].i;
 330     }
 331
 332     for (int j=0;j<3;j++) {
 333       double sh[]=sh_force[j];
 334       for (int i=0;i<mdsize;i++) {
 335         sh[i] = sh[i] * hsq2;
 336       }
 337     }
 338   }
 339
 340   public void run() {
 341     /* Parameter determination */
 342
 343     int mdsize;
 344     double tmpden;
 345     int movemx=50;
 346     //Barrier barr=new Barrier("128.195.175.84");
 347     particle[] one;
 348     int id;
 349     id=this.id;
 350     mdsize = mymd.PARTSIZE;
 351     one=new particle[mdsize];
 352     l = mymd.LENGTH;
 353     tmpden = den;
 354     side = Math.pow((mdsize/tmpden),0.3333333);
 355     rcoff = mm/4.0;
 356
 357     a = side/mm;
 358     sideh = side*0.5;
 359     hsq = h*h;
 360     hsq2 = hsq*0.5;
 361     npartm = mdsize - 1;
 362     rcoffs = rcoff * rcoff;
 363     tscale = 16.0 / (1.0 * mdsize - 1.0);
 364     vaver = 1.13 * Math.sqrt(tref / 24.0);
 365     vaverh = vaver * h;
 366
 367     /* Particle Generation */
 368
 369     xvelocity = 0.0;
 370     yvelocity = 0.0;
 371     zvelocity = 0.0;
 372     ijk = 0;
 373     init(one, mdsize);
 374
 375
 376     /* Synchronise threads and start timer before MD simulation */
 377
 378     //Barrier.enterBarrier(barr);
 379
 380     /* MD simulation */
 381
 382     for (int move=0;move<movemx;move++) {
 383       /* move the particles and update velocities */
 384       for (int i=0;i<mdsize;i++) {
 385         one[i].domove(side,i);
 386       }
 387
 388       /* Barrier */
 389       //Barrier.enterBarrier(barr);
 390
 391       if(id==0) {
 392         doinit(mdsize);
 393       }
 394
 395       mymd.epot[id].d = 0.0;
 396       mymd.vir[id].d = 0.0;
 397       mymd.interacts[id].i = 0;
 398
 399
 400       /* Barrier */
 401       //Barrier.enterBarrier(barr);
 402
 403       /* compute forces */
 404
 405       for (int i=0+id;i<mdsize;i+=nthreads) {
 406         one[i].force(side,rcoff,mdsize,i,xx,yy,zz,mymd);
 407       }
 408
 409       /* Barrier */
 410       //Barrier.enterBarrier(barr);
 411
 412       /* update force arrays */
 413       if(id == 0) {
 414         doinit2(mdsize);
 415       }
 416
 417       /* Barrier */
 418       //Barrier.enterBarrier(barr);
 419
 420       /*scale forces, update velocities */
 421       sum = 0.0;
 422       for (int i=0;i<mdsize;i++) {
 423         sum = sum + one[i].mkekin(hsq2,i);
 424       }
 425
 426       ekin = sum/hsq;
 427
 428       vel = 0.0;
 429       count = 0.0;
 430
 431       /* average velocity */
 432
 433       for (int i=0;i<mdsize;i++) {
 434         velt = one[i].velavg(vaverh,h);
 435         if(velt > vaverh) { count = count + 1.0; }
 436         vel = vel + velt;
 437       }
 438
 439       vel = vel / h;
 440
 441       /* temperature scale if required */
 442
 443       if((move < istop) && (((move+1) % irep) == 0)) {
 444         sc = Math.sqrt(tref / (tscale*ekin));
 445         for (int i=0;i<mdsize;i++) {
 446           one[i].dscal(sc,1);
 447         }
 448         ekin = tref / tscale;
 449       }
 450
 451       /* sum to get full potential energy and virial */
 452
 453       if(((move+1) % iprint) == 0) {
 454         mymd.ek[id].d = 24.0*ekin;
 455         mymd.epot[id].d = 4.0*mymd.epot[id].d;
 456         etot = mymd.ek[id].d + mymd.epot[id].d;
 457         temp = tscale * ekin;
 458         pres = tmpden * 16.0 * (ekin - mymd.vir[id].d) / mdsize;
 459         vel = vel / mdsize;
 460         rp = (count / mdsize) * 100.0;
 461       }
 462       //Barrier.enterBarrier(barr);
 463
 464       //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
 465     }
 466   }
 467 }
 468
 469   class particle {
 470
 471     public double xcoord, ycoord, zcoord;
 472     public double xvelocity,yvelocity,zvelocity;
 473     int part_id;
 474     int id;
 475     double [][] sh_force;
 476     double [][][] sh_force2;
 477     particle[] one;
 478
 479     public particle(double xcoord, double ycoord, double zcoord, double xvelocity,
 480         double yvelocity,double zvelocity, double [][] sh_force,
 481         double [][][] sh_force2,int id, particle[] one) {
 482
 483       this.xcoord = xcoord;
 484       this.ycoord = ycoord;
 485       this.zcoord = zcoord;
 486       this.xvelocity = xvelocity;
 487       this.yvelocity = yvelocity;
 488       this.zvelocity = zvelocity;
 489       this.sh_force = sh_force;
 490       this.sh_force2 = sh_force2;
 491       this.id=id;
 492       this.one=one;
 493     }
 494
 495     public void domove(double side,int part_id) {
 496
 497       xcoord = xcoord + xvelocity + sh_force[0][part_id];
 498       ycoord = ycoord + yvelocity + sh_force[1][part_id];
 499       zcoord = zcoord + zvelocity + sh_force[2][part_id];
 500
 501       if(xcoord < 0) { xcoord = xcoord + side; }
 502       if(xcoord > side) { xcoord = xcoord - side; }
 503       if(ycoord < 0) { ycoord = ycoord + side; }
 504       if(ycoord > side) { ycoord = ycoord - side; }
 505       if(zcoord < 0) { zcoord = zcoord + side; }
 506       if(zcoord > side) { zcoord = zcoord - side; }
 507
 508       xvelocity = xvelocity + sh_force[0][part_id];
 509       yvelocity = yvelocity + sh_force[1][part_id];
 510       zvelocity = zvelocity + sh_force[2][part_id];
 511
 512     }
 513
 514     public void force(double side, double rcoff,int mdsize,int x, double xx, double yy, double zz, JGFMolDynBench mymd) {
 515
 516       double sideh;
 517       double rcoffs;
 518
 519       double fxi,fyi,fzi;
 520       double rd,rrd,rrd2,rrd3,rrd4,rrd6,rrd7,r148;
 521       double forcex,forcey,forcez;
 522       int id=this.id;
 523       sideh = 0.5*side;
 524       rcoffs = rcoff*rcoff;
 525
 526       fxi = 0.0;
 527       fyi = 0.0;
 528       fzi = 0.0;
 529
 530       for (int i=x+1;i<mdsize;i++) {
 531         xx = this.xcoord - one[i].xcoord;
 532         yy = this.ycoord - one[i].ycoord;
 533         zz = this.zcoord - one[i].zcoord;
 534
 535         if(xx < (-sideh)) { xx = xx + side; }
 536         if(xx > (sideh))  { xx = xx - side; }
 537         if(yy < (-sideh)) { yy = yy + side; }
 538         if(yy > (sideh))  { yy = yy - side; }
 539         if(zz < (-sideh)) { zz = zz + side; }
 540         if(zz > (sideh))  { zz = zz - side; }
 541
 542
 543         rd = xx*xx + yy*yy + zz*zz;
 544
 545         if(rd <= rcoffs) {
 546           rrd = 1.0/rd;
 547           rrd2 = rrd*rrd;
 548           rrd3 = rrd2*rrd;
 549           rrd4 = rrd2*rrd2;
 550           rrd6 = rrd2*rrd4;
 551           rrd7 = rrd6*rrd;
 552           mymd.epot[id].d += (rrd6 - rrd3);
 553           r148 = rrd7 - 0.5*rrd4;
 554           mymd.vir[id].d += - rd*r148;
 555           forcex = xx * r148;
 556           fxi = fxi + forcex;
 557
 558           sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
 559
 560           forcey = yy * r148;
 561           fyi = fyi + forcey;
 562
 563           sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
 564
 565           forcez = zz * r148;
 566           fzi = fzi + forcez;
 567
 568           sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
 569
 570           mymd.interacts[id].i++;
 571         }
 572
 573       }
 574
 575       sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
 576       sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
 577       sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
 578
 579     }
 580
 581     public double mkekin(double hsq2,int part_id) {
 582
 583       double sumt = 0.0;
 584
 585       xvelocity = xvelocity + sh_force[0][part_id];
 586       yvelocity = yvelocity + sh_force[1][part_id];
 587       zvelocity = zvelocity + sh_force[2][part_id];
 588
 589       sumt = (xvelocity*xvelocity)+(yvelocity*yvelocity)+(zvelocity*zvelocity);
 590       return sumt;
 591     }
 592
 593     public double velavg(double vaverh,double h) {
 594
 595       double velt;
 596       double sq;
 597
 598       sq = Math.sqrt(xvelocity*xvelocity + yvelocity*yvelocity +
 599           zvelocity*zvelocity);
 600
 601       velt = sq;
 602       return velt;
 603     }
 604
 605     public void dscal(double sc,int incx) {
 606       xvelocity = xvelocity * sc;
 607       yvelocity = yvelocity * sc;
 608       zvelocity = zvelocity * sc;
 609     }
 610   }
 611
 612   class random {
 613
 614     public int iseed;
 615     public double v1,v2;
 616
 617     public random(int iseed,double v1,double v2) {
 618       this.iseed = iseed;
 619       this.v1 = v1;
 620       this.v2 = v2;
 621     }
 622
 623     public double update() {
 624
 625       double rand;
 626       double scale= 4.656612875e-10;
 627
 628       int is1,is2,iss2;
 629       int imult=16807;
 630       int imod = 2147483647;
 631
 632       if (iseed<=0) { iseed = 1; }
 633
 634       is2 = iseed % 32768;
 635       is1 = (iseed-is2)/32768;
 636       iss2 = is2 * imult;
 637       is2 = iss2 % 32768;
 638       is1 = (is1*imult+(iss2-is2)/32768) % (65536);
 639
 640       iseed = (is1*32768+is2) % imod;
 641
 642       rand = scale * iseed;
 643
 644       return rand;
 645
 646     }
 647
 648     public double seed() {
 649
 650       double s,u1,u2,r;
 651       s = 1.0;
 652       do {
 653         u1 = update();
 654         u2 = update();
 655
 656         v1 = 2.0 * u1 - 1.0;
 657         v2 = 2.0 * u2 - 1.0;
 658         s = v1*v1 + v2*v2;
 659
 660       } while (s >= 1.0);
 661
 662       r = Math.sqrt(-2.0*Math.log(s)/s);
 663
 664       return r;
 665
 666     }
 667   }
 668
 669