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