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