∇-Nabla: Numerical Analysis BAsed LAnguage

∇-Nabla: Numerical Analysis BAsed LAnguage

nablaLulesh-512.png nablaCoMD-512.png nablaDDFV-512.png nablaMNLDDFV-512.png nablaSethi-512.png nablaSPH-512.png

∇ CoMD

This partial ∇ port comes from an implementation of the Los Alamos National Security benchmark: CoMD, A Molecular Dynamics Proxy Applications Suite.

with particles;

// ****************************************************************************
// * DEFINES
// ****************************************************************************
#define amu_to_m_e (1822.83)
#define bohr_per_atu_to_A_per_s (0.529/2.418e-17)

// ****************************************************************************
// * OPTIONS
// ****************************************************************************
options{
   option_nx     = 20;
   option_ny     = 20;
   option_nz     = 20;
   option_lat       = 2.07514e-317;
   option_defgrad   = 1.0; // stretch in x direction (deformation gradient)
   option_cutoff    = 4.59;
   option_boxfactor = 1.0;
   option_max_iteration = 20;
};

// ****************************************************************************
// * Déclaration des variables aux mailles (alias 'boxes')
// ****************************************************************************
cells{
   natom;
  ℝ³ dcenter;
   cell_energy;
   cell_stress;
};

// ****************************************************************************
// * Déclaration des variables aux particules
// ****************************************************************************
particles{
   id;
   iType;  // the type of the atom
  ℝ³ r;        // position
  ℝ³ μ;        // momenta of the atom
  ℝ³ f;
   energy;
};

// ****************************************************************************
// * Déclaration des variables globales
// ****************************************************************************
global{
   xNbx;   // number of boxes in each dimension
   yNbx;   // number of boxes in each dimension
   zNbx;   // number of boxes in each dimension
   nboxes; // total number of boxes
   ntot;   // total number of atoms
   xBoxsize;  // size of domains
   yBoxsize;  // size of domains
   zBoxsize;  // size of domains
   rho;       // rhosum for EAM potential
   fi;        // rhobar for EAM potential
   e;         // the total energy of the system
   stress;    // virial stress
   defgrad;   // deformation gradient
   lat;       // Lattice constant
   xBounds;   // x periodic bounds
   yBounds;   // y periodic bounds
   zBounds;   // z periodic bounds
   boxfactor;
  // Lennard-Jones potential (pmd_base_potential_t *pot)
   σ;       // the finite distance at which the inter-particle potential is zero
   ε;       // the depth of the potential well
   cutoff;  // potential cutoff distance in Angstroms 
   mass;    // mass of the atom in atomic mass units
  // Simulation needs
   ts;
   te;
   nsteps;
   ns;
   s6;
   etot;
   r2cut;
};

// ********************************************************
// * fonctions outils
// ********************************************************
 usPeratom( tend,  tstart){
  return ()(()((1.0e6*(tend-tstart)/()ns/()ntot)*1000.0))/1000.0;
}

 timeNow(void){
  return ()(()(ElapsedTime*1000.0))/1000.0;
}

// ****************************************************************************
// * Partie d'initialisation ]-∞,-0[
// ****************************************************************************

// ********************************************************
// * récupération des options
// ********************************************************
void iniOptions(void) @ -12.0{
  cutoff    = option_cutoff;
  defgrad   = option_defgrad;
  boxfactor = option_boxfactor;
}

// ********************************************************
// * Initialisation propre au potentiel Lennard-Jones
// ********************************************************
void iniLJPotentialValues(void) @ -11.0{
  σ=1.53;
  ε=0.0085;
  cutoff=3.0*σ;
  mass = 1.0;
  lat=1.122462048*cutoff;
  info()<<"[iniLJPotentialValues] cutoff="<<cutoff<<", lattice="<<lat;
}

// ********************************************************
// * Initialisation propre au potentiel Lennard-Jones
// ********************************************************
void iniLattice(void) @ -10.0{
  xBounds = option_nx * lat * defgrad;
  yBounds = option_ny * lat;
  zBounds = option_nz * lat;
}

// ********************************************************
// * Initialisation du domaine
// ********************************************************
void iniDomain(void) @ -9.0{
   strain =1.0;
  info() << "bounds = "<<xBounds<<","<<yBounds<<","<<zBounds
         << ", cutoff = "<<cutoff<<", box factor = "<<boxfactor
         << ", strain = "<<strain;
  xNbx = ()floor(xBounds/(cutoff*boxfactor*strain));
  yNbx = ()floor(yBounds/(cutoff*boxfactor*strain));
  zNbx = ()floor(zBounds/(cutoff*boxfactor*strain));
  nboxes = xNbx*yNbx*zNbx;
  xBoxsize=xBounds/()xNbx;
  yBoxsize=yBounds/()yNbx;
  zBoxsize=zBounds/()zNbx;
}

// ********************************************************
// * Initialisation des particules
// ********************************************************
 particles void iniParticles(void) out (particle id) @ -9.0{
  f=μ=r=0.0;
  energy=0.0;
  id=iType=0;
}

// ********************************************************
// * Initialisation des mailles et calcul du centre
// ********************************************************
 cells void computeDomainCenter(void) @ -8.0{
  natom=0.0;
  dcenter=0.0;
   node dcenter+=coord;
  dcenter/=nbNode;
}

// given x,y,z in world co-ordinates return the box in which those coordinates fall
 getBoxIDWorldCoords( x,  y,  z){
   ibox,ibx[3];
  ibx[0] = ()(floor(x/xBoxsize));
  ibx[1] = ()(floor(y/yBoxsize));
  ibx[2] = ()(floor(z/zBoxsize));
  ibox = ibx[0]+xNbx*ibx[1]+xNbx*yNbx*ibx[2];
  return ibox;
}

// finds an appropriate box for an atom based on the spatial coordinates and puts it in there
void putAtomInBox(const  n,  x,  y,  z) {
  ℝ³ ρ;
  // push atom into primary period
  if (x<0.0) x+=xBounds; else if (x>=xBounds) x-=xBounds;
  if (y<0.0) y+=yBounds; else if (y>=yBounds) y-=yBounds;
  if (z<0.0) z+=zBounds; else if (z>=zBounds) z-=zBounds;
  ρ.x=x; ρ.y=y; ρ.z=z;
  particleAddToCell(n,getBoxIDWorldCoords(x,y,z),ρ);
  ntot=ntot+1;
}

// ********************************************************
// * Remplissage des boîtes
// ********************************************************
void putAtomsInBoxes(void) @ -7.0{
   x,y,z;
   i,j,k,n;
   halflat=lat/2.0;
  for(z=lat/4.0,i=j=k=n=0; z<zBounds; z+=halflat,k++)
    for(y=lat/4.0; y<yBounds; y+=halflat,j++)
      for(x=lat*defgrad/4.0; x<xBounds; x+=halflat*defgrad,i++)
        if ((i+j+k) % 2) putAtomInBox(n++,x,y,z);
  particleSyncToCell();
}


// ********************************************************
// * Initialisation des particules
// ********************************************************
 particles void iniAtomTheBox(void) in (cell dcenter)
  out (particle id, particle iType, particle r) @ -6.0{
  id = uid;
  iType = 1;
  r -= dcenter;
}


// ********************************************************
// * Mise à jour du nombre de particules dans chaque maille
// ********************************************************
 own cells void iniParticlesInCells(void) in (particle id, particle r)@ -5.9{
   particle{
    natom+=1.0;
    //cell_energy=0.0;
    //if (!(uid%1024)) info()<<"cell #"<<uid<<" has particle "<< id <<" @ coords="<<r;
  }
}


// ********************************************************
// * Initialisation des variables de la simulation
// ********************************************************
void initSim(void) @ -5.0{
  info() << "total atoms is: "<< ntot;
  info() << "box factor is "<< xBoxsize/cutoff
         << " " << yBoxsize/cutoff << " " << zBoxsize/cutoff;
  nsteps=10;
  info() << "Starting simulation";
  δt = 1.0e-15;
  ts=timeNow();
  ns=()(1+nsteps);
}


// ****************************************************************************
// * calculates forces for the 12-6 lennard jones potential
// * Notes on LJ:
// * http://en.wikipedia.org/wiki/Lennard_Jones_potential
// * LJ is a simple potential of the form:
// * e_lj(r) = 4*ε*((σ/r)^12 - (σ/r)^6)
// * F(r) = 4*ε*(12*σ^12/r^13 - 6*σ^6/r^7)
// * ε and σ are the adjustable parameters in the potential.
// *    ε = well depth
// *    σ   = hard sphere diameter
// * You can also adjust the 12 & 6, but few people do.
// * Some typical values for ε and σ:
// *   material            ε            σ   
// *     N                 36.20 K             3.299 A
// *     O                 44.06 K             2.956 A
// ****************************************************************************
void iniComputeForce(void){
  e = 0.0;
  etot = 0.0;
  stress = 0.0;
  s6 = σ*σ*σ*σ*σ*σ;
  r2cut = cutoff*cutoff;
}
 particles void particleZeroForcesEnergy(void){
  f=0.0;
  energy=0.0;
}
void cellComputeForce(Cell i, Cell j){
   etotl=0.0;
   strss=0.0;
   r2,r6,rs,fr;
  ℝ³ dr;
  ℝ³ drbox=dcenter[i]-dcenter[j];
   i particle{
    strss -= - μ[pi].x*μ[pi].x/mass;
     j particle{
      if (id[pj]<=id[pi]) continue;
      dr = drbox+(r[pi]-r[pj]);
      r2 = dot(dr,dr);
      if (r2>r2cut) continue;
      // Important note: from this point on r actually refers to 1.0/r
      r2 = 1.0/r2;
      r6 = (r2*r2*r2);
      rs = 0.5*r6*(s6*r6 - 1.0);
      energy[pi] += rs; 
      energy[pj] += rs;
      etotl += r6*(s6*r6 - 1.0);
      // different formulation to avoid sqrt computation
      fr = 4.0*ε*s6*r6*r2*(12.0*s6*r6-6.0);
      f[pi] += dr*fr;
      f[pj] -= dr*fr;
      strss += 2.0*fr*dr.x*dr.x;
    } // loop over atoms in j
  } // loop over atoms in i
  cell_energy[i]+=etotl;
  cell_stress[i]+=strss;
}
// loop over all boxes in system via pairs
 own cells void cellsComputeForceViaPairs(void) {
  cell_stress=cell_energy=0.0;
  cellComputeForce(*this,*this);            // itself
  foreach cell cellComputeForce(*this,*cc);  // neighbors
}
 own cells  computeSumEnergy( nrj) in (cell cell_energy){
  nrj+=cell_energy;
  return nrj;
}
own cells  computeSumStress( s) in (cell cell_stress){
  s+=cell_stress;
  return s;
}
void endComputeForce(void) {
  etot = etot*4.0*ε*s6;
  e = etot;
  // renormalize stress
  stress = stress/(xBounds*yBounds*zBounds);
}
void computeForce(void){
  iniComputeForce();
  particleZeroForcesEnergy();
  cellsComputeForceViaPairs();
  etot=computeSumEnergy(etot);
  stress=computeSumStress(stress);
  etot=mpi_reduce(ReduceSum,etot);
  stress=mpi_reduce(ReduceSum,stress);
  endComputeForce();
}

// ********************************************************
// * Initialisation des forces avant la boucle
// ********************************************************
void firstComputeForce(void) @ -0.5{
  computeForce();
}

// ****************************************************************************
// * Partie de calcul ]+0,+oo[ 
// ****************************************************************************

// ****************************************************************************
// * do_compute_work
// ****************************************************************************
void do_compute_work(void) @ 1.0{
  te=timeNow();
  info() << "Iteration #"<<(GlobalIteration-1)*nsteps<<", total system energy e="<< e
         << ", computed in "<<(te-ts)<<"s ("<< usPeratom(te,ts) <<"us/atom for "<<ntot<<" atoms)";
  ts=timeNow();
}

// ****************************************************************************
// * advanceVelocity
// ****************************************************************************
 particles void advanceVelocity( dt) in (particle f) out (particle p){
  μ -= dt*f;
}

// ****************************************************************************
// * All computations are done in atomic units.
// * (http://en.wikipedia.org/wiki/Atomic_units)
// * In these units the mass of an electrion is
// * assumed to be 1.0.  Since our mass is in
// * amu, we need to multiply it by 1822.83 to
// * get the mass in amu.
// ****************************************************************************
 particles void advancePositions( dt) in (particle mass){
  r += dt*μ/(amu_to_m_e*mass);
}

 getIx( iBox) { return (iBox%xNbx);}
 getIy( iBox) { return (iBox/xNbx)%(yNbx);}
 getIz( iBox) { return (iBox/xNbx/yNbx)%zNbx;}
 getIBoxFromIxyz3( x,  y,  z){
   rtn;
  info()<<"\33[7m getIBoxFromIxyz3 "<<x<<"x"<<y<<"x"<<z<<"\33[m";
  if (x<0) x+=xNbx; if (x>=xNbx) x-=xNbx;
  if (y<0) y+=yNbx; if (y>=yNbx) y-=yNbx;
  if (z<0) z+=zNbx; if (z>=zNbx) z-=zNbx;
  info()<<"\33[7m getIBoxFromIxyz3 "<<x<<"x"<<y<<"x"<<z<<"\33[m";
  rtn = x + y*xNbx + z*xNbx*yNbx;
  info()<<"\33[7m getIBoxFromIxyz3 rtn="<<rtn<<"\33[m";
  return rtn;
}

// ****************************************************************************
// * reBoxAll
// ****************************************************************************
 cells void reBoxAll(void){
   ibox=uid;
   iXold=getIx(uid);
   iYold=getIy(uid);
   iZold=getIz(uid);
  //info()<<"\33[7mreBoxAll cell #"<<uid<<":"<<iXold<<"x"<<iYold<<"x"<<iZold<<"\33[m";
  foreach particle{
     jbox;
     iXnew,iYnew,iZnew;
    ℝ³ rnew=r;

    info()<<"\33[35m In box #"<<uid<<", particle #"<<id<<" r="<<r<<"\33[m";
    info()<<" xBoxsize="<<xBoxsize()<<", yBoxsize="<<yBoxsize()<<", zBoxsize="<<zBoxsize();

    if (r.x < 0.0) { iXnew = iXold-1; rnew.x += xBoxsize; }
    else if (r.x >= xBoxsize) { iXnew = iXold+1; rnew.x -= xBoxsize; }
    else { iXnew = iXold; }

    if (r.y < 0.0) { iYnew = iYold-1; rnew.y += yBoxsize; }
    else if (r.y >= yBoxsize) { iYnew = iYold+1; rnew.y -= yBoxsize; }
    else { iYnew = iYold; }

    if (r.z < 0.0) { iZnew = iZold-1; rnew.z += zBoxsize; }
    else if (r.z >= zBoxsize) { iZnew = iZold+1; rnew.z -= zBoxsize; }
    else { iZnew = iZold; }

    jbox = getIBoxFromIxyz3(iXnew,iYnew,iZnew);

    if((jbox<0)||(jbox==ibox)){
      /* do nothing if same box or non-periodic boundary */
    }else{
      r=rnew;      
      moveAtomInBox(id, jbox);
    }
  }
}

// ****************************************************************************
// * Standard verlet algorithm:
// *   1: advance positions half a timestep using current velocities
// *   2: compute forces
// *   3: advance velocities (momenta) a full timestep
// *   4: advance positions half a timestep to bring in sync with velocities.
// ****************************************************************************
void nTimeSteps(void) @ 2.0{
   i;
   dt = δt * bohr_per_atu_to_A_per_s;
  for(i=0;i<nsteps;i++){
    advancePositions(dt/2.);
    computeForce();
    advanceVelocity(dt); 
    advancePositions(dt/2.);    
  }
  // compute force to make consistent
  computeForce();
  te=timeNow();
}

// ****************************************************************************
// * exit_compute_work
// ****************************************************************************
void exit_compute_work(void) @ 3.0{
  if (GlobalIteration > option_max_iteration) exit;
}

camierjs@nabla-lang.org, Jean-Sylvain.Camier@cea.fr