Commit e76cad35 authored by Tim Mattox's avatar Tim Mattox
Browse files

Merge branch 'FixRX_Kokkos_merge' into USER-DPD_kokkos

parents 0a751c59 acc5bde0

src/KOKKOS/Install.sh

+2 −0
Original line number Diff line number Diff line
@@ -109,6 +109,8 @@ action fix_wall_reflect_kokkos.cpp 
action fix_wall_reflect_kokkos.h

action fix_dpd_energy_kokkos.cpp fix_dpd_energy.cpp

action fix_dpd_energy_kokkos.h fix_dpd_energy.h

action fix_rx_kokkos.cpp fix_rx.cpp

action fix_rx_kokkos.h fix_rx.h

action gridcomm_kokkos.cpp gridcomm.cpp

action gridcomm_kokkos.h gridcomm.h

action improper_class2_kokkos.cpp improper_class2.cpp

src/KOKKOS/fix_rx_kokkos.cpp

0 → 100644
+2262 −0

File added.

Preview size limit exceeded, changes collapsed.

src/KOKKOS/fix_rx_kokkos.h

0 → 100644
+280 −0
Original line number Diff line number Diff line 
/* ----------------------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov



   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.



   See the README file in the top-level LAMMPS directory.

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



#ifdef FIX_CLASS



FixStyle(rx/kk,FixRxKokkos<LMPDeviceType>)

FixStyle(rx/kk/device,FixRxKokkos<LMPDeviceType>)

FixStyle(rx/kk/host,FixRxKokkos<LMPHostType>)



#else



#ifndef LMP_FIX_RX_KOKKOS_H

#define LMP_FIX_RX_KOKKOS_H



#include "fix_rx.h"

#include "pair_dpd_fdt_energy_kokkos.h"

#include "kokkos_type.h"

#include "neigh_list.h"

#include "neigh_list_kokkos.h"



namespace LAMMPS_NS {



struct Tag_FixRxKokkos_zeroTemperatureViews {};

struct Tag_FixRxKokkos_zeroCounterViews {};



template <int WT_FLAG, bool NEWTON_PAIR, int NEIGHFLAG>

struct Tag_FixRxKokkos_firstPairOperator {};



template <int WT_FLAG, int LOCAL_TEMP_FLAG>

struct Tag_FixRxKokkos_2ndPairOperator {};



template <bool ZERO_RATES>

struct Tag_FixRxKokkos_solveSystems {};



struct s_CounterType

{

  int nSteps, nIters, nFuncs, nFails;



  KOKKOS_INLINE_FUNCTION

  s_CounterType() : nSteps(0), nIters(0), nFuncs(0), nFails(0) {};



  KOKKOS_INLINE_FUNCTION

  s_CounterType& operator+=(const s_CounterType &rhs)

  {

    nSteps += rhs.nSteps;

    nIters += rhs.nIters;

    nFuncs += rhs.nFuncs;

    nFails += rhs.nFails;

    return *this;

  }



  KOKKOS_INLINE_FUNCTION

  volatile s_CounterType& operator+=(const volatile s_CounterType &rhs) volatile

  {

    nSteps += rhs.nSteps;

    nIters += rhs.nIters;

    nFuncs += rhs.nFuncs;

    nFails += rhs.nFails;

    return *this;

  }

};

typedef struct s_CounterType CounterType;



template <typename DeviceType>

class FixRxKokkos : public FixRX {

 public:

  FixRxKokkos(class LAMMPS *, int, char **);

  virtual ~FixRxKokkos();

  virtual void init();

  void init_list(int, class NeighList *);

  void post_constructor();

  virtual void setup_pre_force(int);

  virtual void pre_force(int);



  // Define a value_type here for the reduction operator on CounterType.

  typedef CounterType value_type;



  KOKKOS_INLINE_FUNCTION

  void operator()(Tag_FixRxKokkos_zeroCounterViews, const int&) const;



  KOKKOS_INLINE_FUNCTION

  void operator()(Tag_FixRxKokkos_zeroTemperatureViews, const int&) const;



  template <int WT_FLAG, bool NEWTON_PAIR, int NEIGHFLAG>

  KOKKOS_INLINE_FUNCTION

  void operator()(Tag_FixRxKokkos_firstPairOperator<WT_FLAG,NEWTON_PAIR,NEIGHFLAG>, const int&) const;



  template <int WT_FLAG, int LOCAL_TEMP_FLAG>

  KOKKOS_INLINE_FUNCTION

  void operator()(Tag_FixRxKokkos_2ndPairOperator<WT_FLAG,LOCAL_TEMP_FLAG>, const int&) const;



  template <bool ZERO_RATES>

  KOKKOS_INLINE_FUNCTION

  void operator()(Tag_FixRxKokkos_solveSystems<ZERO_RATES>, const int&, CounterType&) const;



 //protected:

  PairDPDfdtEnergyKokkos<DeviceType>* pairDPDEKK;

  double VDPD;



  double boltz;

  double t_stop;



  template <typename T, int stride = 1>

  struct StridedArrayType

  {

    typedef T value_type;

    enum { Stride = stride };



    value_type *m_data;



    KOKKOS_INLINE_FUNCTION

    StridedArrayType() : m_data(NULL) {}

    KOKKOS_INLINE_FUNCTION

    StridedArrayType(value_type *ptr) : m_data(ptr) {}



    KOKKOS_INLINE_FUNCTION       value_type& operator()(const int idx)       { return m_data[Stride*idx]; }

    KOKKOS_INLINE_FUNCTION const value_type& operator()(const int idx) const { return m_data[Stride*idx]; }

    KOKKOS_INLINE_FUNCTION       value_type& operator[](const int idx)       { return m_data[Stride*idx]; }

    KOKKOS_INLINE_FUNCTION const value_type& operator[](const int idx) const { return m_data[Stride*idx]; }

  };



  template <int stride = 1>

  struct UserRHSDataKokkos

  {

    StridedArrayType<double,1> kFor;

    StridedArrayType<double,1> rxnRateLaw;

  };



  void solve_reactions(const int vflag, const bool isPreForce);



  int rhs       (double, const double *, double *, void *) const;

  int rhs_dense (double, const double *, double *, void *) const;

  int rhs_sparse(double, const double *, double *, void *) const;



  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  int k_rhs       (double, const VectorType&, VectorType&, UserDataType& ) const;



  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  int k_rhs_dense (double, const VectorType&, VectorType&, UserDataType& ) const;



  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  int k_rhs_sparse(double, const VectorType&, VectorType&, UserDataType& ) const;



  //!< Classic Runge-Kutta 4th-order stepper.

  void rk4(const double t_stop, double *y, double *rwork, void *v_params) const;



  //!< Runge-Kutta-Fehlberg ODE Solver.

  void rkf45(const int neq, const double t_stop, double *y, double *rwork, void *v_params, CounterType& counter) const;



  //!< Runge-Kutta-Fehlberg ODE stepper function.

  void rkf45_step (const int neq, const double h, double y[], double y_out[],

                   double rwk[], void *) const;



  //!< Initial step size estimation for the Runge-Kutta-Fehlberg ODE solver.

  int rkf45_h0 (const int neq, const double t, const double t_stop,

                     const double hmin, const double hmax,

                     double& h0, double y[], double rwk[], void *v_params) const;



  //!< Classic Runge-Kutta 4th-order stepper.

  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  void k_rk4(const double t_stop, VectorType& y, VectorType& rwork, UserDataType& userData) const;



  //!< Runge-Kutta-Fehlberg ODE Solver.

  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  void k_rkf45(const int neq, const double t_stop, VectorType& y, VectorType& rwork, UserDataType& userData, CounterType& counter) const;



  //!< Runge-Kutta-Fehlberg ODE stepper function.

  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  void k_rkf45_step (const int neq, const double h, VectorType& y, VectorType& y_out,

                     VectorType& rwk, UserDataType& userData) const;



  //!< Initial step size estimation for the Runge-Kutta-Fehlberg ODE solver.

  template <typename VectorType, typename UserDataType>

    KOKKOS_INLINE_FUNCTION

  int k_rkf45_h0 (const int neq, const double t, const double t_stop,

                  const double hmin, const double hmax,

                  double& h0, VectorType& y, VectorType& rwk, UserDataType& userData) const;



  //!< ODE Solver diagnostics.

  void odeDiagnostics(void);



  //!< Special counters per-ode.

  int *diagnosticCounterPerODEnSteps;

  int *diagnosticCounterPerODEnFuncs;

  DAT::tdual_int_1d k_diagnosticCounterPerODEnSteps;

  DAT::tdual_int_1d k_diagnosticCounterPerODEnFuncs;

  //typename ArrayTypes<DeviceType>::t_int_1d d_diagnosticCounterPerODEnSteps;

  //typename ArrayTypes<DeviceType>::t_int_1d d_diagnosticCounterPerODEnFuncs;

  typename DAT::t_int_1d d_diagnosticCounterPerODEnSteps;

  typename DAT::t_int_1d d_diagnosticCounterPerODEnFuncs;

  typename HAT::t_int_1d h_diagnosticCounterPerODEnSteps;

  typename HAT::t_int_1d h_diagnosticCounterPerODEnFuncs;



  template <typename KokkosDeviceType>

  struct KineticsType

  {

    // Arrhenius rate coefficients.

    typename ArrayTypes<KokkosDeviceType>::t_float_1d Arr, nArr, Ea;



    // Dense versions.

    typename ArrayTypes<KokkosDeviceType>::t_float_2d stoich, stoichReactants, stoichProducts;



    // Sparse versions.

    typename ArrayTypes<KokkosDeviceType>::t_int_2d   nuk, inu;

    typename ArrayTypes<KokkosDeviceType>::t_float_2d nu;

    typename ArrayTypes<KokkosDeviceType>::t_int_1d   isIntegral;

  };



  //!< Kokkos versions of the kinetics data.

  KineticsType<LMPHostType> h_kineticsData;

  KineticsType<DeviceType>  d_kineticsData;



  bool update_kinetics_data;



  void create_kinetics_data(void);



  // Need a dual-view and device-view for dpdThetaLocal and sumWeights since they're used in several callbacks.

  DAT::tdual_efloat_1d k_dpdThetaLocal, k_sumWeights;

  //typename ArrayTypes<DeviceType>::t_efloat_1d d_dpdThetaLocal, d_sumWeights;

  typename DAT::t_efloat_1d d_dpdThetaLocal, d_sumWeights;

  typename HAT::t_efloat_1d h_dpdThetaLocal, h_sumWeights;



  typename ArrayTypes<DeviceType>::t_x_array_randomread d_x       ;

  typename ArrayTypes<DeviceType>::t_int_1d_randomread  d_type    ;

  typename ArrayTypes<DeviceType>::t_efloat_1d          d_dpdTheta;



  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;

  typename ArrayTypes<DeviceType>::t_ffloat_2d     d_cutsq;

  //double **h_cutsq;



  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;

  typename ArrayTypes<DeviceType>::t_int_1d       d_ilist    ;

  typename ArrayTypes<DeviceType>::t_int_1d       d_numneigh ;



  typename ArrayTypes<DeviceType>::t_float_2d  d_dvector;

  typename ArrayTypes<DeviceType>::t_int_1d    d_mask   ;



  typename ArrayTypes<DeviceType>::t_double_1d d_scratchSpace;

  size_t scratchSpaceSize;



  // Error flag for any failures.

  DAT::tdual_int_scalar k_error_flag;



  template <int WT_FLAG, int LOCAL_TEMP_FLAG, bool NEWTON_PAIR, int NEIGHFLAG>

  void computeLocalTemperature();



  int pack_reverse_comm(int, int, double *);

  void unpack_reverse_comm(int, int *, double *);

  int pack_forward_comm(int , int *, double *, int, int *);

  void unpack_forward_comm(int , int , double *);



 //private: // replicate a few from FixRX

  int my_restartFlag;

  int nlocal;

};



}



#endif

#endif



/* ERROR/WARNING messages:



*/

src/USER-DPD/fix_rx.cpp

+179 −68
Original line number Diff line number Diff line
@@ -220,6 +220,9 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) : 


FixRX::~FixRX()

{

  //printf("Inside FixRX::~FixRX copymode= %d\n", copymode);

  if (copymode) return;



  // De-Allocate memory to prevent memory leak

  for (int ii = 0; ii < nreactions; ii++){

    delete [] stoich[ii];
@@ -673,7 +676,17 @@ void FixRX::setup_pre_force(int vflag) 


  if(restartFlag){

    restartFlag = 0;

  } else {

  }

  else

  {

    int ode_counter[4] = {0};



    UserRHSData userData;

    userData.kFor = new double[nreactions];

    userData.rxnRateLaw = new double[nreactions];



    double *rwork = new double[8*nspecies];



    if(localTempFlag){

      int count = nlocal + (newton_pair ? nghost : 0);

      dpdThetaLocal = new double[count];
@@ -686,22 +699,27 @@ void FixRX::setup_pre_force(int vflag) 
        tmp = atom->dvector[ispecies][id];

        atom->dvector[ispecies+nspecies][id] = tmp;

      }



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

      if (mask[i] & groupbit){



        // Set the reaction rate constants to zero:  no reactions occur at step 0

        for(int irxn=0;irxn<nreactions;irxn++)

          kR[irxn] = 0.0;

          userData.kFor[irxn] = 0.0;



        if (odeIntegrationFlag == ODE_LAMMPS_RK4)

          rk4(i,NULL);

          rk4(i, rwork, &userData);

        else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)

          rkf45(i,NULL);

          rkf45(i, rwork, &userData, ode_counter);

      }



    // Communicate the updated momenta and velocities to all nodes

    comm->forward_comm_fix(this);

    if(localTempFlag) delete [] dpdThetaLocal;



    delete [] userData.kFor;

    delete [] userData.rxnRateLaw;

    delete [] rwork;

  }

}
@@ -709,12 +727,13 @@ void FixRX::setup_pre_force(int vflag) 


void FixRX::pre_force(int vflag)

{

  TimerType timer_start = getTimeStamp();



  int nlocal = atom->nlocal;

  int nghost = atom->nghost;

  int *mask = atom->mask;

  double *dpdTheta = atom->dpdTheta;

  int newton_pair = force->newton_pair;

  double theta;



  if(localTempFlag){

    int count = nlocal + (newton_pair ? nghost : 0);
@@ -726,7 +745,10 @@ void FixRX::pre_force(int vflag) 
  TimerType timer_localTemperature = getTimeStamp();



  // Zero the counters for the ODE solvers.

  this->nSteps = this->nIters = this->nFuncs = this->nFails = 0;

  int nSteps = 0;

  int nIters = 0;

  int nFuncs = 0;

  int nFails = 0;



  if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency == 1)

  {
@@ -734,10 +756,23 @@ void FixRX::pre_force(int vflag) 
    memory->create( diagnosticCounterPerODE[FuncSum], nlocal, "FixRX::diagnosticCounterPerODE");

  }



  double *rwork = new double[8*nspecies + nreactions];

  //#pragma omp parallel \

  //   reduction(+: nSteps, nIters, nFuncs, nFails )

  {

    double *rwork = new double[8*nspecies];



    UserRHSData userData;

    userData.kFor = new double[nreactions];

    userData.rxnRateLaw = new double[nreactions];



    int ode_counter[4] = { 0 };



    //#pragma omp for schedule(runtime)

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

    if (mask[i] & groupbit){

    {

      if (mask[i] & groupbit)

      {

        double theta;

        if (localTempFlag)

          theta = dpdThetaLocal[i];

        else
@@ -745,24 +780,42 @@ void FixRX::pre_force(int vflag) 


        //Compute the reaction rate constants

        for (int irxn = 0; irxn < nreactions; irxn++)

        kR[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);

          userData.kFor[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);



        if (odeIntegrationFlag == ODE_LAMMPS_RK4)

        rk4(i,rwork);

          rk4(i, rwork, &userData);

        else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)

        rkf45(i,rwork);

          rkf45(i, rwork, &userData, ode_counter);

      }

    }



  TimerType timer_ODE = getTimeStamp();

    nSteps += ode_counter[0];

    nIters += ode_counter[1];

    nFuncs += ode_counter[2];

    nFails += ode_counter[3];



    delete [] rwork;

    delete [] userData.kFor;

    delete [] userData.rxnRateLaw;



  } // end parallel region



  TimerType timer_ODE = getTimeStamp();



  // Communicate the updated momenta and velocities to all nodes

  comm->forward_comm_fix(this);

  if(localTempFlag) delete [] dpdThetaLocal;



  TimerType timer_stop = getTimeStamp();



  double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE);



  //printf("me= %d total= %g temp= %g ode= %g comm= %g nlocal= %d nfc= %d %d\n", comm->me,

  //                       getElapsedTime(timer_start, timer_stop),

  //                       getElapsedTime(timer_start, timer_localTemperature),

  //                       getElapsedTime(timer_localTemperature, timer_ODE),

  //                       getElapsedTime(timer_ODE, timer_stop), nlocal, nFuncs, nSteps);



  // Warn the user if a failure was detected in the ODE solver.

  if (nFails > 0){

    char sbuf[128];
@@ -958,21 +1011,15 @@ void FixRX::setupParams() 


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



void FixRX::rk4(int id, double *rwork)

void FixRX::rk4(int id, double *rwork, void* v_params)

{

  double *k1 = NULL;

  if (rwork == NULL)

    k1 = new double[6*nspecies + nreactions];

  else

    k1 = rwork;

  double *k1 = rwork;

  double *k2 = k1 + nspecies;

  double *k3 = k2 + nspecies;

  double *k4 = k3 + nspecies;

  double *y  = k4 + nspecies;

  double *yp = y  + nspecies;



  double *dummyArray = yp + nspecies; // Passed to the rhs function.



  const int numSteps = minSteps;



  const double h = update->dt / double(numSteps);
@@ -989,25 +1036,25 @@ void FixRX::rk4(int id, double *rwork) 
  for (int step = 0; step < numSteps; step++)

  {

    // k1

    rhs(0.0,y,k1,dummyArray);

    rhs(0.0,y,k1,v_params);



    // k2

    for (int ispecies = 0; ispecies < nspecies; ispecies++)

      yp[ispecies] = y[ispecies] + 0.5*h*k1[ispecies];



    rhs(0.0,yp,k2,dummyArray);

    rhs(0.0,yp,k2,v_params);



    // k3

    for (int ispecies = 0; ispecies < nspecies; ispecies++)

      yp[ispecies] = y[ispecies] + 0.5*h*k2[ispecies];



    rhs(0.0,yp,k3,dummyArray);

    rhs(0.0,yp,k3,v_params);



    // k4

    for (int ispecies = 0; ispecies < nspecies; ispecies++)

      yp[ispecies] = y[ispecies] + h*k3[ispecies];



    rhs(0.0,yp,k4,dummyArray);

    rhs(0.0,yp,k4,v_params);



    for (int ispecies = 0; ispecies < nspecies; ispecies++)

      y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 + k3[ispecies]/3.0 + k4[ispecies]/6.0);
@@ -1022,9 +1069,6 @@ void FixRX::rk4(int id, double *rwork) 
      y[ispecies] = 0.0;

    atom->dvector[ispecies][id] = y[ispecies];

  }



  if (rwork == NULL)

    delete [] k1;

}



/* ---------------------------------------------------------------------- */
@@ -1274,6 +1318,78 @@ void FixRX::odeDiagnostics(void) 
  double max_per_proc[numCounters];

  double min_per_proc[numCounters];



  if(1)

  {

     static bool firstStep = true;



     static TimerType oldTimeStamp (-1);



     TimerType now = getTimeStamp();



     // Query the fix database and look for rx_weight for the balance fix.

     int type_flag = -1;

     int rx_weight_index = atom->find_custom( "rx_weight", /*0:int, 1:float*/ type_flag );



     // Compute the average # of neighbors.

     double averageNumNeighbors = 0;

     {

        const int inum = pairDPDE->list->inum;

        const int* ilist = pairDPDE->list->ilist;

        const int* numneigh = pairDPDE->list->numneigh;



        for (int ii = 0; ii < inum; ++ii)

        {

           const int i = ilist[ii];

           averageNumNeighbors += numneigh[i];

        }



        averageNumNeighbors /= inum;

     }



     printf("me= %d nst= %g nfc= %g time= %g nlocal= %g lmpnst= %g weight_idx= %d 1st= %d aveNeigh= %g\n", comm->me, this->diagnosticCounter[0], this->diagnosticCounter[1], this->diagnosticCounter[2], this->diagnosticCounter[3], this->diagnosticCounter[4], rx_weight_index, firstStep, averageNumNeighbors);



     if (rx_weight_index != -1 && !firstStep && 0)

     {

        double *rx_weight = atom->dvector[rx_weight_index];



        const int nlocal = atom->nlocal;

        const int *mask = atom->mask;



        if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency == 1)

        {

          const double total_time = getElapsedTime( oldTimeStamp, now );

          const double fixrx_time = this->diagnosticCounter[TimeSum];

          const double time_ratio = fixrx_time / total_time;



          double tsum = 0.0;

          double tmin = 100000, tmax = 0;

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

            if (mask[i] & groupbit)

            {

              double nfunc_ratio = double( diagnosticCounterPerODE[FuncSum][i] ) / diagnosticCounter[FuncSum];

              rx_weight[i] = nfunc_ratio * fixrx_time + (total_time - fixrx_time) / nlocal;

              tmin = fmin( tmin, rx_weight[i] );

              tmax = fmax( tmax, rx_weight[i] );

              tsum += rx_weight[i];

              //rx_weight[i] = (double) diagnosticCounterPerODE[FuncSum][i];

            }



          printf("me= %d total= %g fixrx= %g ratio= %g tsum= %g %g %g %g\n", comm->me, total_time, fixrx_time, time_ratio, tsum, (total_time - fixrx_time) / nlocal, tmin, tmax);

        }

        else

        {

          error->warning(FLERR, "Dynamic load balancing enabled but per-atom weights not available.");



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

            if (mask[i] & groupbit)

              rx_weight[i] = 1.0;

        }

     }



     firstStep = false;

     oldTimeStamp = now;

  }



  // Compute counters per dpd time-step.

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

    my_vals[i] = this->diagnosticCounter[i] / nTimes;
@@ -1347,7 +1463,7 @@ void FixRX::odeDiagnostics(void) 
    if (screen)  fprintf(screen,"%s\n", smesg); \

    if (logfile) fprintf(logfile,"%s\n", smesg); }



    sprintf(smesg, "FixRX::ODE Diagnostics:  # of steps  |# of rhs evals| run-time (sec)");

    sprintf(smesg, "FixRX::ODE Diagnostics:  # of iters  |# of rhs evals| run-time (sec) | # atoms");

    print_mesg(smesg);



    sprintf(smesg, "         AVG per ODE  : %-12.5g | %-12.5g | %-12.5g", avg_per_atom[0], avg_per_atom[1], avg_per_atom[2]);
@@ -1369,7 +1485,7 @@ void FixRX::odeDiagnostics(void) 
      print_mesg(smesg);

    }



    sprintf(smesg, "         AVG per Proc : %-12.5g | %-12.5g | %-12.5g", avg_per_proc[0], avg_per_proc[1], avg_per_proc[2]);

    sprintf(smesg, "         AVG per Proc : %-12.5g | %-12.5g | %-12.5g | %-12.5g", avg_per_proc[StepSum], avg_per_proc[FuncSum], avg_per_proc[TimeSum], avg_per_proc[AtomSum]);

    print_mesg(smesg);



    if (comm->nprocs > 1){
@@ -1377,13 +1493,13 @@ void FixRX::odeDiagnostics(void) 
      for (int i = 0; i < numCounters; ++i)

        rms_per_proc[i] = sqrt( sum_sq[i] / comm->nprocs );



      sprintf(smesg, "         RMS per Proc : %-12.5g | %-12.5g | %-12.5g", rms_per_proc[0], rms_per_proc[1], rms_per_proc[2]);

      sprintf(smesg, "         RMS per Proc : %-12.5g | %-12.5g | %-12.5g | %-12.5g", rms_per_proc[0], rms_per_proc[1], rms_per_proc[2], rms_per_proc[AtomSum]);

      print_mesg(smesg);



      sprintf(smesg, "         MAX per Proc : %-12.5g | %-12.5g | %-12.5g", max_per_proc[0], max_per_proc[1], max_per_proc[2]);

      sprintf(smesg, "         MAX per Proc : %-12.5g | %-12.5g | %-12.5g | %-12.5g", max_per_proc[0], max_per_proc[1], max_per_proc[2], max_per_proc[AtomSum]);

      print_mesg(smesg);



      sprintf(smesg, "         MIN per Proc : %-12.5g | %-12.5g | %-12.5g", min_per_proc[0], min_per_proc[1], min_per_proc[2]);

      sprintf(smesg, "         MIN per Proc : %-12.5g | %-12.5g | %-12.5g | %-12.5g", min_per_proc[0], min_per_proc[1], min_per_proc[2], min_per_proc[AtomSum]);

      print_mesg(smesg);

    }
@@ -1403,7 +1519,7 @@ void FixRX::odeDiagnostics(void) 
  return;

}



void FixRX::rkf45(int id, double *rwork)

void FixRX::rkf45(int id, double *rwork, void *v_param, int ode_counter[])

{

  // Rounding coefficient.

  const double uround = DBL_EPSILON;
@@ -1412,12 +1528,7 @@ void FixRX::rkf45(int id, double *rwork) 
  const double adaption_limit = 4.0;



  //double *y = new double[8*nspecies + nreactions];

  double *y = NULL;

  if (rwork == NULL)

    y = new double[8*nspecies + nreactions];

  else

    y = rwork;

  double *rhstmp = y + 8*nspecies;

  double *y = rwork;



  const int neq = nspecies;
@@ -1454,7 +1565,7 @@ void FixRX::rkf45(int id, double *rwork) 
  if (h < h_min){

    //fprintf(stderr,"hin not implemented yet\n");

    //exit(-1);

    nfe = rkf45_h0 (neq, t, t_stop, h_min, h_max, h, y, y + neq, rhstmp);

    nfe = rkf45_h0 (neq, t, t_stop, h_min, h_max, h, y, y + neq, v_param);

  }



  //printf("t= %e t_stop= %e h= %e\n", t, t_stop, h);
@@ -1465,7 +1576,7 @@ void FixRX::rkf45(int id, double *rwork) 
    double *eout = yout + neq;



    // Take a trial step.

    rkf45_step (neq, h, y, yout, eout, rhstmp);

    rkf45_step (neq, h, y, yout, eout, v_param);



    // Estimate the solution error.

      // ... weighted 2-norm of the error.
@@ -1513,16 +1624,17 @@ void FixRX::rkf45(int id, double *rwork) 


    if (maxIters && nit > maxIters){

      //fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop);

      nFails ++;

      //nFails ++;

      ode_counter[3] ++;

      break;

      // We should set an error here so that the solution is not used!

    }



  } // end while



  nSteps += nst;

  nIters += nit;

  nFuncs += nfe;

  ode_counter[0] += nst;

  ode_counter[1] += nit;

  ode_counter[2] += nfe;



  //if (diagnosticFrequency == 1 && diagnosticCounterPerODE[StepSum] != NULL)

  if (diagnosticCounterPerODE[StepSum] != NULL){
@@ -1539,9 +1651,6 @@ void FixRX::rkf45(int id, double *rwork) 
      y[ispecies] = 0.0;

    atom->dvector[ispecies][id] = y[ispecies];

  }



  if (rwork == NULL)

    delete [] y;

}



/* ---------------------------------------------------------------------- */
@@ -1559,21 +1668,23 @@ int FixRX::rhs(double t, const double *y, double *dydt, void *params) 


int FixRX::rhs_dense(double t, const double *y, double *dydt, void *params)

{

  double rxnRateLawForward;

  double *rxnRateLaw = (double *) params;

  double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;

  double concentration;

  int nspecies = atom->nspecies_dpd;

  UserRHSData *userData = (UserRHSData *) params;



  double *rxnRateLaw = userData->rxnRateLaw;

  double *kFor       = userData->kFor;



  const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;

  const int nspecies = atom->nspecies_dpd;



  for(int ispecies=0; ispecies<nspecies; ispecies++)

    dydt[ispecies] = 0.0;



  // Construct the reaction rate laws

  for(int jrxn=0; jrxn<nreactions; jrxn++){

    rxnRateLawForward = kR[jrxn];

    double rxnRateLawForward = kFor[jrxn];



    for(int ispecies=0; ispecies<nspecies; ispecies++){

      concentration = y[ispecies]/VDPD;

      const double concentration = y[ispecies]/VDPD;

      rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);

    }

    rxnRateLaw[jrxn] = rxnRateLawForward;
@@ -1591,13 +1702,13 @@ int FixRX::rhs_dense(double t, const double *y, double *dydt, void *params) 


int FixRX::rhs_sparse(double t, const double *y, double *dydt, void *v_params) const

{

   double *_rxnRateLaw = (double *) v_params;

   UserRHSData *userData = (UserRHSData *) v_params;



   const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;



   #define kFor         (this->kR)

   #define kFor         (userData->kFor)

   #define kRev         (NULL)

   #define rxnRateLaw   (_rxnRateLaw)

   #define rxnRateLaw   (userData->rxnRateLaw)

   #define conc         (dydt)

   #define maxReactants (this->sparseKinetics_maxReactants)

   #define maxSpecies   (this->sparseKinetics_maxSpecies)

src/USER-DPD/fix_rx.h

+15 −8

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE