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