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[4];
  73     mid[0] = (128<<24)|(195<<16)|(175<<8)|79;
  74     mid[1] = (128<<24)|(195<<16)|(175<<8)|80;
  75     mid[2] = (128<<24)|(195<<16)|(175<<8)|78;
  76     mid[3] = (128<<24)|(195<<16)|(175<<8)|73;
  77
  78     double sh_force [][];
  79     double sh_force2 [][][];
  80     int partsize, numthreads;
  81     atomic {
  82       partsize = mold.PARTSIZE;
  83       numthreads = mold.nthreads;
  84       mybarr = global new BarrierServer(numthreads);
  85     }
  86     mybarr.start(mid[0]);
  87
  88     atomic {
  89       sh_force = global new double[3][partsize];
  90       sh_force2 = global new double[3][numthreads][partsize];
  91       mold.epot = global new DoubleWrapper[numthreads];
  92       mold.vir  = global new DoubleWrapper[numthreads];
  93       mold.ek   = global new DoubleWrapper[numthreads];
  94       mold.interacts = global new IntWrapper[numthreads];
  95       for(int i=0;i<numthreads;i++) {
  96         mold.epot[i]=global new DoubleWrapper();
  97         mold.vir[i]=global new DoubleWrapper();
  98         mold.ek[i]=global new DoubleWrapper();
  99         mold.interacts[i]=global new IntWrapper();
 100       }
 101     }
 102
 103     // spawn threads
 104     MDWrap[] thobjects = new MDWrap[numthreads];
 105
 106     atomic {
 107       for(int i=0;i<numthreads;i++) {
 108         thobjects[i] = new MDWrap(global new mdRunner(i,mold.mm,sh_force,sh_force2,mold.nthreads,mold));
 109       }
 110     }
 111
 112     boolean waitfordone=true;
 113     while(waitfordone) {
 114       atomic{
 115         if (mybarr.done)
 116           waitfordone=false;
 117       }
 118     }
 119
 120     for(int i=0;i<numthreads;i++) {
 121       thobjects[i].md.start(mid[i]);
 122     }
 123
 124     for(int i=0;i<numthreads;i++) {
 125       thobjects[i].md.join();
 126     }
 127   }
 128
 129   public void JGFvalidate(){
 130     double[] refval = new double[2];
 131     refval[0] = 1731.4306625334357;
 132     refval[1] = 7397.392307839352;
 133     double dev = Math.fabs(ek[0].d - refval[size]);
 134     if (dev > 1.0e-10 ){
 135       //System.printString("Validation failed\n");
 136       //System.printString("Kinetic Energy = " + (long)ek[0] + "  " + (long)dev + "  " + size + "\n");
 137     }
 138   }
 139 }
 140
 141 class mdRunner extends Thread {
 142
 143   double count;
 144   int id,i,j,k,lg,mm;
 145   double l,rcoff,rcoffs,side,sideh,hsq,hsq2,vel,velt;
 146   double a,r,sum,tscale,sc,ekin,ts,sp;
 147   double den;
 148   double tref;
 149   double h;
 150   double vaver,vaverh,rand;
 151   double etot,temp,pres,rp;
 152   double u1, u2, s, xx, yy, zz;
 153   double xvelocity, yvelocity, zvelocity;
 154
 155   double [][] sh_force;
 156   double [][][] sh_force2;
 157
 158   int ijk,npartm,iseed,tint;
 159   int irep;
 160   int istop;
 161   int iprint;
 162
 163   JGFMolDynBench mymd;
 164   int nthreads;
 165
 166   public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,
 167                   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.nthreads = nthreads;
 173     this.mymd = mymd;
 174     count = 0.0;
 175     den = 0.83134;
 176     tref = 0.722;
 177     h = 0.064;
 178     irep = 10;
 179     istop = 19;
 180     iprint = 10;
 181   }
 182
 183     public void init(particle[] one, int mdsize) {
 184         int id=this.id;
 185         for (int lg=0; lg<=1; lg++) {
 186         for (int i=0; i<mm; i++) {
 187           for (int j=0; j<mm; j++) {
 188             for (int k=0; k<mm; k++) {
 189               one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
 190                   xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
 191               ijk = ijk + 1;
 192             }
 193           }
 194         }
 195       }
 196
 197       for (int lg=1; lg<=2; lg++) {
 198         for (int i=0; i<mm; i++) {
 199           for (int j=0; j<mm; j++) {
 200             for (int k=0; k<mm; k++) {
 201               one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
 202                                       (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
 203               ijk = ijk + 1;
 204             }
 205           }
 206         }
 207       }
 208
 209       /* Initialise velocities */
 210
 211       iseed = 0;
 212       double v1 = 0.0;
 213       double v2 = 0.0;
 214       random randnum = new random(iseed,v1,v2);
 215
 216
 217       for (int i=0; i<mdsize; i+=2) {
 218         r  = randnum.seed();
 219         one[i].xvelocity = r*randnum.v1;
 220         one[i+1].xvelocity  = r*randnum.v2;
 221       }
 222
 223       for (int i=0; i<mdsize; i+=2) {
 224         r  = randnum.seed();
 225         one[i].yvelocity = r*randnum.v1;
 226         one[i+1].yvelocity  = r*randnum.v2;
 227       }
 228
 229       for (int i=0; i<mdsize; i+=2) {
 230         r  = randnum.seed();
 231         one[i].zvelocity = r*randnum.v1;
 232         one[i+1].zvelocity  = r*randnum.v2;
 233       }
 234
 235
 236       /* velocity scaling */
 237
 238       ekin = 0.0;
 239       sp = 0.0;
 240
 241       for(int i=0;i<mdsize;i++) {
 242         sp = sp + one[i].xvelocity;
 243       }
 244       sp = sp / mdsize;
 245
 246       for(int i=0;i<mdsize;i++) {
 247         one[i].xvelocity = one[i].xvelocity - sp;
 248         ekin = ekin + one[i].xvelocity*one[i].xvelocity;
 249       }
 250
 251       sp = 0.0;
 252       for(int i=0;i<mdsize;i++) {
 253         sp = sp + one[i].yvelocity;
 254       }
 255       sp = sp / mdsize;
 256
 257       for(int i=0;i<mdsize;i++) {
 258         one[i].yvelocity = one[i].yvelocity - sp;
 259         ekin = ekin + one[i].yvelocity*one[i].yvelocity;
 260       }
 261
 262
 263       sp = 0.0;
 264       for(int i=0;i<mdsize;i++) {
 265         sp = sp + one[i].zvelocity;
 266       }
 267       sp = sp / mdsize;
 268
 269       for(int i=0;i<mdsize;i++) {
 270         one[i].zvelocity = one[i].zvelocity - sp;
 271         ekin = ekin + one[i].zvelocity*one[i].zvelocity;
 272       }
 273
 274       ts = tscale * ekin;
 275       sc = h * Math.sqrt(tref/ts);
 276
 277
 278       for(int i=0;i<mdsize;i++) {
 279
 280         one[i].xvelocity = one[i].xvelocity * sc;
 281         one[i].yvelocity = one[i].yvelocity * sc;
 282         one[i].zvelocity = one[i].zvelocity * sc;
 283
 284       }
 285
 286     }
 287
 288     public void doinit(int mdsize) {
 289         for(int j=0;j<3;j++) {
 290             double[] sh=sh_force[j];
 291             for (int i=0;i<mdsize;i++) {
 292                 sh[i] = 0.0;
 293             }
 294         }
 295     }
 296
 297
 298     public void doinit2(int mdsize) {
 299           for(int k=0;k<3;k++) {
 300               double[] sh=sh_force[k];
 301               double [][] sha=sh_force2[k];
 302               for(int j=0;j<nthreads;j++) {
 303                   double[] sha2=sha[j];
 304                   for(int i=0;i<mdsize;i++) {
 305                 sh[i] += sha2[i];
 306               }
 307             }
 308           }
 309
 310           for(int k=0;k<3;k++) {
 311               double [][] sh1=sh_force2[k];
 312               for(int j=0;j<nthreads;j++) {
 313                   double[] sh2=sh1[j];
 314             for(int i=0;i<mdsize;i++) {
 315
 316
 317                 sh2[i] = 0.0;
 318               }
 319             }
 320           }
 321
 322           for(int j=1;j<nthreads;j++) {
 323             mymd.epot[0].d += mymd.epot[j].d;
 324             mymd.vir[0].d += mymd.vir[j].d;
 325           }
 326           for(int j=1;j<nthreads;j++) {
 327             mymd.epot[j].d = mymd.epot[0].d;
 328             mymd.vir[j].d = mymd.vir[0].d;
 329           }
 330           for(int j=0;j<nthreads;j++) {
 331             mymd.interactions += mymd.interacts[j].i;
 332           }
 333
 334           for (int j=0;j<3;j++) {
 335               double sh[]=sh_force[j];
 336             for (int i=0;i<mdsize;i++) {
 337               sh[i] = sh[i] * hsq2;
 338             }
 339           }
 340 }
 341
 342   public void run() {
 343     /* Parameter determination */
 344
 345     int mdsize;
 346     double tmpden;
 347     int movemx=50;
 348     Barrier barr=new Barrier("128.195.175.79");
 349     particle[] one;
 350     int id;
 351     atomic {
 352         id=this.id;
 353       mdsize = mymd.PARTSIZE;
 354       one=new particle[mdsize];
 355       l = mymd.LENGTH;
 356       tmpden = den;
 357       side = Math.pow((mdsize/tmpden),0.3333333);
 358       rcoff = mm/4.0;
 359
 360       a = side/mm;
 361       sideh = side*0.5;
 362       hsq = h*h;
 363       hsq2 = hsq*0.5;
 364       npartm = mdsize - 1;
 365       rcoffs = rcoff * rcoff;
 366       tscale = 16.0 / (1.0 * mdsize - 1.0);
 367       vaver = 1.13 * Math.sqrt(tref / 24.0);
 368       vaverh = vaver * h;
 369
 370       /* Particle Generation */
 371
 372       xvelocity = 0.0;
 373       yvelocity = 0.0;
 374       zvelocity = 0.0;
 375       ijk = 0;
 376       init(one, mdsize);
 377
 378     }
 379
 380     /* Synchronise threads and start timer before MD simulation */
 381
 382     Barrier.enterBarrier(barr);
 383
 384     /* MD simulation */
 385
 386     for (int move=0;move<movemx;move++) {
 387       atomic {
 388         /* move the particles and update velocities */
 389         for (int i=0;i<mdsize;i++) {
 390           one[i].domove(side,i);
 391         }
 392       }
 393
 394       /* Barrier */
 395       Barrier.enterBarrier(barr);
 396
 397       atomic {
 398         if(id==0) {
 399             doinit(mdsize);
 400         }
 401
 402         mymd.epot[id].d = 0.0;
 403         mymd.vir[id].d = 0.0;
 404         mymd.interacts[id].i = 0;
 405       }
 406
 407
 408       /* Barrier */
 409       Barrier.enterBarrier(barr);
 410
 411       atomic {
 412         /* compute forces */
 413
 414         for (int i=0+id;i<mdsize;i+=nthreads) {
 415           one[i].force(side,rcoff,mdsize,i,xx,yy,zz,mymd);
 416         }
 417       }
 418
 419       /* Barrier */
 420       Barrier.enterBarrier(barr);
 421
 422       /* update force arrays */
 423       atomic {
 424         if(id == 0) {
 425             doinit2(mdsize);
 426         }
 427       }
 428
 429       /* Barrier */
 430       Barrier.enterBarrier(barr);
 431
 432       atomic {
 433         /*scale forces, update velocities */
 434         sum = 0.0;
 435         for (int i=0;i<mdsize;i++) {
 436           sum = sum + one[i].mkekin(hsq2,i);
 437         }
 438
 439         ekin = sum/hsq;
 440
 441         vel = 0.0;
 442         count = 0.0;
 443
 444         /* average velocity */
 445
 446         for (int i=0;i<mdsize;i++) {
 447           velt = one[i].velavg(vaverh,h);
 448           if(velt > vaverh) { count = count + 1.0; }
 449           vel = vel + velt;
 450         }
 451
 452         vel = vel / h;
 453
 454         /* temperature scale if required */
 455
 456         if((move < istop) && (((move+1) % irep) == 0)) {
 457           sc = Math.sqrt(tref / (tscale*ekin));
 458           for (int i=0;i<mdsize;i++) {
 459             one[i].dscal(sc,1);
 460           }
 461           ekin = tref / tscale;
 462         }
 463
 464         /* sum to get full potential energy and virial */
 465
 466         if(((move+1) % iprint) == 0) {
 467           mymd.ek[id].d = 24.0*ekin;
 468           mymd.epot[id].d = 4.0*mymd.epot[id].d;
 469           etot = mymd.ek[id].d + mymd.epot[id].d;
 470           temp = tscale * ekin;
 471           pres = tmpden * 16.0 * (ekin - mymd.vir[id].d) / mdsize;
 472           vel = vel / mdsize;
 473           rp = (count / mdsize) * 100.0;
 474         }
 475       }
 476       Barrier.enterBarrier(barr);
 477     }
 478
 479     //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
 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   global double [][] sh_force;
 491   global double [][][] sh_force2;
 492   particle[] one;
 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, particle[] one) {
 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.one=one;
 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     int id=this.id;
 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 - one[i].xcoord;
 547       yy = this.ycoord - one[i].ycoord;
 548       zz = this.zcoord - 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].d += (rrd6 - rrd3);
 568         r148 = rrd7 - 0.5*rrd4;
 569         mymd.vir[id].d += - 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].i++;
 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     xvelocity = xvelocity * sc;
 622     yvelocity = yvelocity * sc;
 623     zvelocity = zvelocity * sc;
 624   }
 625 }
 626
 627 class random {
 628
 629   public int iseed;
 630   public double v1,v2;
 631
 632   public random(int iseed,double v1,double v2) {
 633     this.iseed = iseed;
 634     this.v1 = v1;
 635     this.v2 = v2;
 636   }
 637
 638   public double update() {
 639
 640     double rand;
 641     double scale= 4.656612875e-10;
 642
 643     int is1,is2,iss2;
 644     int imult=16807;
 645     int imod = 2147483647;
 646
 647     if (iseed<=0) { iseed = 1; }
 648
 649     is2 = iseed % 32768;
 650     is1 = (iseed-is2)/32768;
 651     iss2 = is2 * imult;
 652     is2 = iss2 % 32768;
 653     is1 = (is1*imult+(iss2-is2)/32768) % (65536);
 654
 655     iseed = (is1*32768+is2) % imod;
 656
 657     rand = scale * iseed;
 658
 659     return rand;
 660
 661   }
 662
 663   public double seed() {
 664
 665     double s,u1,u2,r;
 666     s = 1.0;
 667     do {
 668       u1 = update();
 669       u2 = update();
 670
 671       v1 = 2.0 * u1 - 1.0;
 672       v2 = 2.0 * u2 - 1.0;
 673       s = v1*v1 + v2*v2;
 674
 675     } while (s >= 1.0);
 676
 677     r = Math.sqrt(-2.0*Math.log(s)/s);
 678
 679     return r;
 680
 681   }
 682 }
 683
 684