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