Commit 605fe53c authored by Tim Mattox's avatar Tim Mattox
Browse files

USER-DPD: make the kokkos and non-kokkos fix_shardlow match more closely 

This was accomplished with several key changes:
1) Modified fix_shardlow's control flow to match fix_shardlow_kokkos so
that random numbers are pulled fromn the RNGs in exactly the same order.

2) Created random_external_state.h, a simplified version of the Kokkos
random number generator that keeps its state variables external to itself.
Thus it can be used both with and without Kokkos enabled, as long as the
caller stores and passes in the required state variable.

3) Replaced all references to random_mars.h and Kokkos_Random.hpp code in
the fix_shardlow* files with calls to the random_external_state.h code,
guaranteeing that fix_shardlow* is using an identical RNG in all cases.

Result: most (56 of 61) of our internal tests now generate the same results
with kokkos turned on or off.  Four cases still differ due to what appear
to be vectorization caused rounding differences, and the fifth case
appears to be something triggered by the kokkos "atom_style hybrid" code.
parent 39df9f5d

src/.gitignore

+1 −0
Original line number Diff line number Diff line
@@ -546,6 +546,7 @@ 
/fix_shake.h

/fix_shardlow.cpp

/fix_shardlow.h

/random_external_state.h

/fix_smd.cpp

/fix_smd.h

/fix_species.cpp

src/KOKKOS/fix_shardlow_kokkos.cpp

+17 −58
Original line number Diff line number Diff line
@@ -61,6 +61,7 @@ 


using namespace LAMMPS_NS;

using namespace FixConst;

using namespace random_external_state;



#define EPSILON 1.0e-10

#define EPSILON_SQUARED ((EPSILON) * (EPSILON))
@@ -89,10 +90,6 @@ FixShardlowKokkos<DeviceType>::FixShardlowKokkos(LAMMPS *lmp, int narg, char **a 
//   if(k_pairDPDE){

    comm_forward = 3;

    comm_reverse = 5;

    maxRNG = 0;

#ifdef DPD_USE_RAN_MARS

    pp_random = NULL;

#endif

//   } else {

//     comm_forward = 3;

//     comm_reverse = 3;
@@ -121,13 +118,6 @@ template<class DeviceType> 
FixShardlowKokkos<DeviceType>::~FixShardlowKokkos()

{

  ghostmax = 0;

#ifdef DPD_USE_RAN_MARS

  if (pp_random) {

    for (int i = 1; i < maxRNG; ++i) delete pp_random[i];

    delete[] pp_random;

    pp_random = NULL;

  }

#endif

}



/* ---------------------------------------------------------------------- */
@@ -279,11 +269,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd( 
  int start_ii, int count, int id

)

{

#ifdef DPD_USE_RAN_MARS

  class RanMars *pRNG = pp_random[id];

#else

  rand_type rand_gen = rand_pool.get_state(id);

#endif

  es_RNG_t RNGstate = d_rand_state(id);



  int ct = count;

  int ii = start_ii;
@@ -345,12 +331,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd( 
        double halfsigma_ij = STACKPARAMS?m_params[itype][jtype].halfsigma:params(itype,jtype).halfsigma;

        double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;



        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v *

#ifdef DPD_USE_RAN_MARS

            pRNG->gaussian();

#else

            rand_gen.normal();

#endif

        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);



        const double mass_j = masses(massPerI ? j : jtype);

        double massinv_j = 1.0 / mass_j;
@@ -412,9 +393,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd( 
    v(i, 2) = vzi;

  }



#ifndef DPD_USE_RAN_MARS

  rand_pool.free_state(rand_gen);

#endif

  d_rand_state(id) = RNGstate;

}

#endif
@@ -431,11 +410,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde( 
  int start_ii, int count, int id

) const

{

#ifdef DPD_USE_RAN_MARS

  class RanMars *pRNG = pp_random[id];

#else

  rand_type rand_gen = rand_pool.get_state(id);

#endif

  es_RNG_t RNGstate = d_rand_state(id);



  int ct = count;

  int ii = start_ii;
@@ -506,12 +481,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde( 
        double halfsigma_ij = STACKPARAMS?m_params[itype][jtype].halfsigma:params(itype,jtype).halfsigma;

        double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;



        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v *

#ifdef DPD_USE_RAN_MARS

            pRNG->gaussian();

#else

            rand_gen.normal();

#endif

        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);



        const double mass_j = masses(massPerI ? j : jtype);

        double mass_ij_div_neg4_ftm2v = mass_j*mass_i_div_neg4_ftm2v;
@@ -520,12 +490,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde( 
        // Compute uCond

        double kappa_ij = STACKPARAMS?m_params[itype][jtype].kappa:params(itype,jtype).kappa;

        double alpha_ij = STACKPARAMS?m_params[itype][jtype].alpha:params(itype,jtype).alpha;

        double del_uCond = alpha_ij*wr*dtsqrt *

#ifdef DPD_USE_RAN_MARS

            pRNG->gaussian();

#else

            rand_gen.normal();

#endif

        double del_uCond = alpha_ij*wr*dtsqrt * es_normal(RNGstate);



        del_uCond += kappa_ij*(theta_i_inv - theta_j_inv)*wdt;

        uCond[j] -= del_uCond;
@@ -601,9 +566,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde( 
    ii++;

  }



#ifndef DPD_USE_RAN_MARS

  rand_pool.free_state(rand_gen);

#endif

  d_rand_state(id) = RNGstate;

}
@@ -644,20 +607,16 @@ void FixShardlowKokkos<DeviceType>::initial_integrate(int vflag) 
    maxWorkItemCt = (int) ssa_gitemLoc.dimension_1();

  }

  if (maxWorkItemCt > maxRNG) {

#ifdef DPD_USE_RAN_MARS

    if (pp_random) {

      for (int i = 1; i < maxRNG; ++i) delete pp_random[i];

      delete[] pp_random;

      pp_random = NULL;

    }

    pp_random = new RanMars*[maxWorkItemCt];

    for (int i = 1; i < maxWorkItemCt; ++i) {

      pp_random[i] = new RanMars(lmp, k_pairDPDE->seed + comm->me + comm->nprocs*i);

    es_RNG_t serial_rand_state;

    es_init(serial_rand_state, pairDPDE->seed + comm->me);



    d_rand_state = es_RNGs_type("Kokkos::fix_shardlow::rand_state",maxWorkItemCt);

    typename es_RNGs_type::HostMirror h_rand_state = create_mirror_view(d_rand_state);

    for (int i = 0; i < maxWorkItemCt; ++i) {

      es_genNextParallelState(serial_rand_state, h_rand_state(i));

    }

    pp_random[0] = k_pairDPDE->random;

#else

    rand_pool.init(k_pairDPDE->seed + comm->me, maxWorkItemCt);

#endif

    deep_copy(d_rand_state,h_rand_state);



    maxRNG = maxWorkItemCt;

  }

src/KOKKOS/fix_shardlow_kokkos.h

+4 −11
Original line number Diff line number Diff line
@@ -93,17 +93,6 @@ class FixShardlowKokkos : public FixShardlow { 
#endif

  PairDPDfdtEnergyKokkos<DeviceType> *k_pairDPDE;



  int maxRNG;

#ifdef DPD_USE_RAN_MARS

  class RanMars **pp_random;

#elif defined(DPD_USE_Random_XorShift1024)

  Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;

  typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::generator_type rand_type;

#else

  Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;

  typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;

#endif



  Kokkos::DualView<params_ssa**,Kokkos::LayoutRight,DeviceType> k_params;

  typename Kokkos::DualView<params_ssa**,

    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
@@ -127,6 +116,10 @@ class FixShardlowKokkos : public FixShardlow { 
  typename AT::t_float_1d_randomread masses;

  typename AT::t_efloat_1d dpdTheta;



  // Storage for the es_RNG state variables

  typedef Kokkos::View<random_external_state::es_RNG_t*,DeviceType> es_RNGs_type;

  es_RNGs_type d_rand_state;



  double dtsqrt; // = sqrt(update->dt);

  int ghostmax;

  int nlocal, nghost;

src/USER-DPD/fix_shardlow.cpp

+76 −40
Original line number Diff line number Diff line
@@ -48,7 +48,6 @@ 
#include "neighbor.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "random_mars.h"

#include "memory.h"

#include "domain.h"

#include "modify.h"
@@ -60,6 +59,7 @@ 


using namespace LAMMPS_NS;

using namespace FixConst;

using namespace random_external_state;



#define EPSILON 1.0e-10

#define EPSILON_SQUARED ((EPSILON) * (EPSILON))
@@ -87,6 +87,7 @@ static const char cite_fix_shardlow[] = 


FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :

  Fix(lmp, narg, arg), pairDPD(NULL), pairDPDE(NULL), v_t0(NULL)

  ,rand_state(NULL)

{

  if (lmp->citeme) lmp->citeme->add(cite_fix_shardlow);
@@ -99,6 +100,7 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) : 
  if (pairDPDE == NULL)

    pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy/kk",1);



  maxRNG = 0;

  if(pairDPDE){

    comm_forward = 3;

    comm_reverse = 5;
@@ -116,6 +118,8 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) : 


FixShardlow::~FixShardlow()

{

  memory->destroy(rand_state);

  maxRNG = 0;

}



/* ---------------------------------------------------------------------- */
@@ -173,12 +177,12 @@ void FixShardlow::setup(int vflag) 
   NOTE: only implemented for orthogonal boxes, not triclinic

------------------------------------------------------------------------- */

void FixShardlow::ssa_update_dpd(

  int i,

  int *jlist,

  int jlen

  int start_ii,

  int count,

  int id

)

{

  class RanMars *pRNG;

  es_RNG_t RNGstate = rand_state[id];

  double **x = atom->x;

  double **v = atom->v;

  double *rmass = atom->rmass;
@@ -192,6 +196,16 @@ void FixShardlow::ssa_update_dpd( 


  const double dt     = update->dt;



  int ct = count;

  int ii = start_ii;



while (ct-- > 0) {

  const int i = list->ilist[ii];

  const int *jlist = list->firstneigh[ii];

  const int jlen = list->numneigh[ii];

  ii++;

  if (jlen <= 0) continue;



  const double xtmp = x[i][0];

  const double ytmp = x[i][1];

  const double ztmp = x[i][2];
@@ -203,7 +217,6 @@ void FixShardlow::ssa_update_dpd( 


  int itype = type[i];



  pRNG = pairDPD->random;

  cut2_i = pairDPD->cutsq[itype];

  cut_i  = pairDPD->cut[itype];

  sigma_i = pairDPD->sigma[itype];
@@ -254,7 +267,7 @@ void FixShardlow::ssa_update_dpd( 
      double halfsigma_ij = 0.5*sigma_i[jtype];

      double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;



      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * pRNG->gaussian();

      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);



      double mass_j = (rmass) ? rmass[j] : mass[jtype];

      double massinv_j = 1.0 / mass_j;
@@ -316,6 +329,9 @@ void FixShardlow::ssa_update_dpd( 
  v[i][2] = vzi;

}



  rand_state[id] = RNGstate;

}



/* ----------------------------------------------------------------------

   Perform the stochastic integration and Shardlow update for constant energy

   Allow for both per-type and per-atom mass
@@ -323,12 +339,12 @@ void FixShardlow::ssa_update_dpd( 
   NOTE: only implemented for orthogonal boxes, not triclinic

------------------------------------------------------------------------- */

void FixShardlow::ssa_update_dpde(

  int i,

  int *jlist,

  int jlen

  int start_ii,

  int count,

  int id

)

{

  class RanMars *pRNG;

  es_RNG_t RNGstate = rand_state[id];

  double **x = atom->x;

  double **v = atom->v;

  double *rmass = atom->rmass;
@@ -346,6 +362,16 @@ void FixShardlow::ssa_update_dpde( 


  const double dt     = update->dt;



  int ct = count;

  int ii = start_ii;



while (ct-- > 0) {

  const int i = list->ilist[ii];

  const int *jlist = list->firstneigh[ii];

  const int jlen = list->numneigh[ii];

  ii++;

  if (jlen <= 0) continue;



  const double xtmp = x[i][0];

  const double ytmp = x[i][1];

  const double ztmp = x[i][2];
@@ -359,7 +385,6 @@ void FixShardlow::ssa_update_dpde( 
  double uCond_i = uCond[i];

  int itype = type[i];



  pRNG = pairDPDE->random;

  cut2_i = pairDPDE->cutsq[itype];

  cut_i  = pairDPDE->cut[itype];

  sigma_i = pairDPDE->sigma[itype];
@@ -415,7 +440,7 @@ void FixShardlow::ssa_update_dpde( 
      double halfsigma_ij = 0.5*sigma_i[jtype];

      double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;



      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * pRNG->gaussian();

      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);



      double mass_j = (rmass) ? rmass[j] : mass[jtype];

      double mass_ij_div_neg4_ftm2v = mass_j*mass_i_div_neg4_ftm2v;
@@ -424,7 +449,7 @@ void FixShardlow::ssa_update_dpde( 
      // Compute uCond

      double kappa_ij = kappa_i[jtype];

      double alpha_ij = sqrt(boltz2*kappa_ij);

      double del_uCond = alpha_ij*wr*dtsqrt * pRNG->gaussian();

      double del_uCond = alpha_ij*wr*dtsqrt * es_normal(RNGstate);



      del_uCond += kappa_ij*(theta_i_inv - theta_j_inv)*wdt;

      uCond[j] -= del_uCond;
@@ -499,6 +524,9 @@ void FixShardlow::ssa_update_dpde( 
  uCond[i] = uCond_i;

}



  rand_state[id] = RNGstate;

}



void FixShardlow::initial_integrate(int vflag)

{

  int i,ii,inum;
@@ -507,7 +535,6 @@ void FixShardlow::initial_integrate(int vflag) 
  int nlocal = atom->nlocal;

  int nghost = atom->nghost;



  int airnum;

  const bool useDPDE = (pairDPDE != NULL);



  // NOTE: this logic is specific to orthogonal boxes, not triclinic
@@ -530,6 +557,30 @@ void FixShardlow::initial_integrate(int vflag) 
    error->one(FLERR, msg);

  }



  NPairHalfBinNewtonSSA *np_ssa = dynamic_cast<NPairHalfBinNewtonSSA*>(list->np);

  if (!np_ssa) error->one(FLERR, "NPair wasn't a NPairHalfBinNewtonSSA object");

  int ssa_phaseCt = np_ssa->ssa_phaseCt;

  int *ssa_phaseLen = np_ssa->ssa_phaseLen;

  int **ssa_itemLoc = np_ssa->ssa_itemLoc;

  int **ssa_itemLen = np_ssa->ssa_itemLen;

  int ssa_gphaseCt = np_ssa->ssa_gphaseCt;

  int *ssa_gphaseLen = np_ssa->ssa_gphaseLen;

  int **ssa_gitemLoc = np_ssa->ssa_gitemLoc;

  int **ssa_gitemLen = np_ssa->ssa_gitemLen;



  int maxWorkItemCt = np_ssa->ssa_maxPhaseLen;

  if (maxWorkItemCt > maxRNG) {

    es_RNG_t serial_rand_state;

    es_init(serial_rand_state, pairDPDE->seed + comm->me);



    memory->grow(rand_state, maxWorkItemCt, "FixShardlow:rand_state");

    for (int i = 0; i < maxWorkItemCt; ++i) {

      es_genNextParallelState(serial_rand_state, rand_state[i]);

    }



    maxRNG = maxWorkItemCt;

  }



#ifdef DEBUG_SSA_PAIR_CT

  for (int i = 0; i < 2; ++i)

    for (int j = 0; j < 3; ++j)
@@ -545,13 +596,6 @@ void FixShardlow::initial_integrate(int vflag) 


  dtsqrt = sqrt(update->dt);



  NPairHalfBinNewtonSSA *np_ssa = dynamic_cast<NPairHalfBinNewtonSSA*>(list->np);

  if (!np_ssa) error->one(FLERR, "NPair wasn't a NPairHalfBinNewtonSSA object");

  int ssa_phaseCt = np_ssa->ssa_phaseCt;

  int *ssa_phaseLen = np_ssa->ssa_phaseLen;

  int **ssa_itemLoc = np_ssa->ssa_itemLoc;

  int **ssa_itemLen = np_ssa->ssa_itemLen;



  // process neighbors in the local AIR

  for (int workPhase = 0; workPhase < ssa_phaseCt; ++workPhase) {

    int workItemCt = ssa_phaseLen[workPhase];
@@ -559,22 +603,14 @@ void FixShardlow::initial_integrate(int vflag) 
    for (int workItem = 0; workItem < workItemCt; ++workItem) {

      int ct = ssa_itemLen[workPhase][workItem];

      ii = ssa_itemLoc[workPhase][workItem];



      while (ct-- > 0) {

        int len = list->numneigh[ii];

        if (len > 0) {

          if (useDPDE) ssa_update_dpde(ilist[ii], list->firstneigh[ii], len);

          else ssa_update_dpd(ilist[ii], list->firstneigh[ii], len);

        }

        ii++;

      }

      if (useDPDE) ssa_update_dpde(ii, ct, workItem);

      else ssa_update_dpd(ii, ct, workItem);

    }

  }



  ii = inum;

  //Loop over all 13 outward directions (7 stages)

  for (airnum = 1; airnum <=7; airnum++){

    int ct = list->AIRct_ssa[airnum];

  for (int workPhase = 0; workPhase < ssa_gphaseCt; ++workPhase) {

    int workItemCt = ssa_gphaseLen[workPhase];



    // Communicate the updated velocities to all nodes

    comm->forward_comm_fix(this);
@@ -585,12 +621,11 @@ void FixShardlow::initial_integrate(int vflag) 
      memset(&(atom->uMech[nlocal]), 0, sizeof(double)*nghost);

    }



    // process neighbors in this AIR

    while (ct-- > 0) {

      int len = list->numneigh[ii];

      if (useDPDE) ssa_update_dpde(ilist[ii], list->firstneigh[ii], len);

      else ssa_update_dpd(ilist[ii], list->firstneigh[ii], len);

      ii++;

    for (int workItem = 0; workItem < workItemCt; ++workItem) {

      int ct = ssa_gitemLen[workPhase][workItem];

      ii = ssa_gitemLoc[workPhase][workItem];

      if (useDPDE) ssa_update_dpde(ii, ct, workItem);

      else ssa_update_dpd(ii, ct, workItem);

    }



    // Communicate the ghost deltas to the atom owners
@@ -699,6 +734,7 @@ double FixShardlow::memory_usage() 
{

  double bytes = 0.0;

  bytes += sizeof(double)*3*atom->nghost; // v_t0[]

  bytes += sizeof(*rand_state)*maxRNG; // rand_state[]

  return bytes;

}

src/USER-DPD/fix_shardlow.h

+6 −2
Original line number Diff line number Diff line
@@ -21,6 +21,8 @@ FixStyle(shardlow,FixShardlow) 
#define LMP_FIX_SHARDLOW_H



#include "fix.h"

#include "random_external_state.h"

#include <math.h>



namespace LAMMPS_NS {
@@ -52,12 +54,14 @@ class FixShardlow : public Fix { 
  class PairDPDfdt *pairDPD;

  class PairDPDfdtEnergy *pairDPDE;

  double (*v_t0)[3];

  int maxRNG;



 private:

  random_external_state::es_RNG_t *rand_state;

  double dtsqrt; // = sqrt(update->dt);



  void ssa_update_dpd(int, int *, int);  // Constant Temperature

  void ssa_update_dpde(int, int *, int); // Constant Energy

  void ssa_update_dpd(int, int, int);  // Constant Temperature

  void ssa_update_dpde(int, int, int); // Constant Energy



};

CRA Git | Maintained and supported by SUSTech CRA and CCSE