Unverified Commit 3bb8294f authored by Steve Plimpton's avatar Steve Plimpton Committed by GitHub
Browse files

Merge pull request #718 from timattox/USER-DPD_es_RNG 

USER-DPD: External State RNG
parents 450c689a 29df5a53

src/.gitignore

+1 −0
Original line number Diff line number Diff line
@@ -547,6 +547,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

+77 −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,31 @@ 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) {

    uint64_t my_seed = comm->me + (useDPDE ? pairDPDE->seed : pairDPD->seed);

    es_RNG_t serial_rand_state;

    es_init(serial_rand_state, my_seed);



    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 +597,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 +604,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 +622,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 +735,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