From cff65b956adcee7cf5c5a6a7868382f9bf17bb41 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 22 Jul 2016 22:52:03 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15355 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-DPD/fix_eos_table_rx.cpp | 2 +- src/USER-DPD/fix_rx.cpp | 32 +++-------------------- src/USER-DPD/fix_rx.h | 2 +- src/USER-DPD/pair_multi_lucy_rx.cpp | 39 +++++++++++------------------ src/dump_image.cpp | 2 +- 5 files changed, 22 insertions(+), 55 deletions(-) diff --git a/src/USER-DPD/fix_eos_table_rx.cpp b/src/USER-DPD/fix_eos_table_rx.cpp index 70199e60ce..1e94c2f87c 100644 --- a/src/USER-DPD/fix_eos_table_rx.cpp +++ b/src/USER-DPD/fix_eos_table_rx.cpp @@ -687,7 +687,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai) // Store the current thetai in t1 t1 = MAX(thetai,tb->lo); - t1 = MIN(thetai,tb->hi); + t1 = MIN(t1,tb->hi); if(t1==tb->hi) delta = -delta; // Compute u1 at thetai diff --git a/src/USER-DPD/fix_rx.cpp b/src/USER-DPD/fix_rx.cpp index 5639646156..96f7f44422 100644 --- a/src/USER-DPD/fix_rx.cpp +++ b/src/USER-DPD/fix_rx.cpp @@ -388,7 +388,7 @@ void FixRX::post_constructor() /* ---------------------------------------------------------------------- */ -int FixRX::initSparse() +void FixRX::initSparse() { const int Verbosity = 1; @@ -637,9 +637,8 @@ void FixRX::setup_pre_force(int vflag) int nlocal = atom->nlocal; int nghost = atom->nghost; int *mask = atom->mask; - double *dpdTheta = atom->dpdTheta; int newton_pair = force->newton_pair; - double tmp, theta; + double tmp; int ii; if(localTempFlag){ @@ -662,8 +661,6 @@ void FixRX::setup_pre_force(int vflag) } for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit){ - if(localTempFlag) theta = dpdThetaLocal[i]; - else theta=dpdTheta[i]; // Set the reaction rate constants to zero: no reactions occur at step 0 for(int irxn=0;irxndpdThetaLocal; } TimerType timer_localTemperature = getTimeStamp(); @@ -746,14 +740,7 @@ void FixRX::pre_force(int vflag) comm->forward_comm_fix(this); if(localTempFlag) delete [] dpdThetaLocal; - TimerType timer_stop = getTimeStamp(); double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE); -// printf("FixRX::pre_force total= %e Temp= %e ODE= %e Comm= %e nlocal= %d %d %d %d %d\n", -// getElapsedTime(timer_start, timer_stop), -// (localTempFlag) ? getElapsedTime(timer_start, timer_localTemperature) : 0.0, -// time_ODE, -// getElapsedTime(timer_ODE, timer_stop), -// nlocal, nSteps, nIters, nFuncs, nFails); // Warn the user if a failure was detected in the ODE solver. if (nFails > 0){ @@ -904,6 +891,7 @@ void FixRX::read_file(char *file) error->all(FLERR,"Illegal fix rx command"); } word = strtok(NULL, " \t\n\r\f"); + if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation"); if(strcmp(word,"=") == 0) sign = 1.0; if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0){ if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation"); @@ -1034,12 +1022,9 @@ void FixRX::rk4(int id, double *rwork) void FixRX::rkf45_step (const int neq, const double h, double y[], double y_out[], double rwk[], void* v_param) { - const double c20=0.25; const double c21=0.25; - const double c30=0.375; const double c31=0.09375; const double c32=0.28125; - const double c40=0.92307692307692; const double c41=0.87938097405553; const double c42=-3.2771961766045; const double c43=3.3208921256258; @@ -1047,19 +1032,16 @@ void FixRX::rkf45_step (const int neq, const double h, double y[], double y_out[ const double c52=-8.0; const double c53=7.1734892787524; const double c54=-0.20589668615984; - const double c60=0.5; const double c61=-0.2962962962963; const double c62=2.0; const double c63=-1.3816764132554; const double c64=0.45297270955166; const double c65=-0.275; const double a1=0.11574074074074; - const double a2=0.0; const double a3=0.54892787524366; const double a4=0.5353313840156; const double a5=-0.2; const double b1=0.11851851851852; - const double b2=0.0; const double b3=0.51898635477583; const double b4=0.50613149034201; const double b5=-0.18; @@ -1146,10 +1128,6 @@ int FixRX::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 double uround = DBL_EPSILON; - const double tdist = fabs(t_stop - t); - const double tround = tdist * uround; - // Set lower and upper bounds on h0, and take geometric mean as first trial value. // Exit with this value if the bounds cross each other. @@ -1696,7 +1674,7 @@ void FixRX::computeLocalTemperature() int newton_pair = force->newton_pair; // local temperature variables - double wij=0.0, fr, fr4; + double wij=0.0; double *dpdTheta = atom->dpdTheta; // Initialize the local density and local temperature arrays @@ -1738,8 +1716,6 @@ void FixRX::computeLocalTemperature() if (rsq < pairDPDE->cutsq[itype][jtype]) { double rcut = sqrt(pairDPDE->cutsq[itype][jtype]); double rij = sqrt(rsq); - double tmpFactor = 1.0-rij/rcut; - double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; double ratio = rij/rcut; // Lucy's Weight Function diff --git a/src/USER-DPD/fix_rx.h b/src/USER-DPD/fix_rx.h index f1b6cb5a4b..fc411e265a 100644 --- a/src/USER-DPD/fix_rx.h +++ b/src/USER-DPD/fix_rx.h @@ -90,7 +90,7 @@ class FixRX : public Fix { // Sparse stoichiometric matrix storage format and methods. bool useSparseKinetics; //SparseKinetics sparseKinetics; - int initSparse(void); + void initSparse(void); int rhs_sparse(double, const double *, double *, void *) const; int sparseKinetics_maxReactants; //type; int nlocal = atom->nlocal; int nghost = atom->nghost; - double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double fractionOld1_i,fractionOld1_j; double fractionOld2_i,fractionOld2_j; - double fraction1_i,fraction1_j; - double fraction2_i,fraction2_j; + double fraction1_i; double *uCG = atom->uCG; double *uCGnew = atom->uCGnew; @@ -166,11 +164,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag) fractionOld1_i = fractionOld1[i]; fractionOld2_i = fractionOld2[i]; fraction1_i = fraction1[i]; - fraction2_i = fraction2[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; @@ -184,8 +180,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag) fractionOld1_j = fractionOld1[j]; fractionOld2_j = fractionOld2[j]; - fraction1_j = fraction1[j]; - fraction2_j = fraction2[j]; tb = &tables[tabindex[itype][jtype]]; if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){ @@ -205,8 +199,10 @@ void PairMultiLucyRX::compute(int eflag, int vflag) } A_i = tb->f[itable]; A_j = tb->f[jtable]; - fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype])); - fpair = fpair/sqrt(rsq); + + const double rfactor = 1.0-sqrt(rsq/cutsq[itype][jtype]); + fpair = 0.5*(A_i + A_j)*(4.0-3.0*rfactor)*rfactor*rfactor*rfactor; + fpair /= sqrt(rsq); } else if (tabstyle == LINEAR) { itable = static_cast ((rho[i]*rho[i] - tb->innersq) * tb->invdelta); @@ -232,9 +228,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag) A_i = tb->f[itable] + fraction_i*tb->df[itable]; A_j = tb->f[jtable] + fraction_j*tb->df[jtable]; - fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype])); - - fpair = fpair / sqrt(rsq); + const double rfactor = 1.0-sqrt(rsq/cutsq[itype][jtype]); + fpair = 0.5*(A_i + A_j)*(4.0-3.0*rfactor)*rfactor*rfactor*rfactor; + fpair /= sqrt(rsq); } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx"); @@ -389,15 +385,11 @@ void PairMultiLucyRX::coeff(int narg, char **arg) // insure cutoff is within table if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length"); - double rlo,rhi; if (tb->rflag == 0) { - rlo = tb->rfile[0]; - rhi = tb->rfile[tb->ninput-1]; + rho_0 = tb->rfile[0]; } else { - rlo = tb->rlo; - rhi = tb->rhi; + rho_0 = tb->rlo; } - rho_0 = rlo; tb->match = 0; if (tabstyle == LINEAR && tb->ninput == tablength && @@ -903,8 +895,8 @@ void PairMultiLucyRX::computeLocalDensity() const double delz = ztmp - x[j][2]; const double rsq = delx*delx + dely*dely + delz*delz; - if (one_type) - if (rsq < cutsq_type11){ + if (one_type) { + if (rsq < cutsq_type11) { const double rcut = rcut_type11; const double r_over_rcut = sqrt(rsq) / rcut; const double tmpFactor = 1.0 - r_over_rcut; @@ -913,9 +905,7 @@ void PairMultiLucyRX::computeLocalDensity() rho_i += factor; if (newton_pair || j < nlocal) rho[j] += factor; - } - else - if (rsq < cutsq[itype][jtype]){ + } else if (rsq < cutsq[itype][jtype]) { const double rcut = sqrt(cutsq[itype][jtype]); const double tmpFactor = 1.0-sqrt(rsq)/rcut; const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; @@ -924,6 +914,7 @@ void PairMultiLucyRX::computeLocalDensity() if (newton_pair || j < nlocal) rho[j] += factor; } + } } rho[i] = rho_i; diff --git a/src/dump_image.cpp b/src/dump_image.cpp index 5d1154478f..fd3dbed9e9 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -1036,7 +1036,7 @@ void DumpImage::create_image() // render objects provided by a fix if (fixflag) { - int tridraw,edgedraw; + int tridraw=0,edgedraw=0; if (domain->dimension == 3) { tridraw = 1; edgedraw = 1;