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