Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80746812
fix_rx.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 2, 09:40
Size
52 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 09:40 (2 d)
Engine
blob
Format
Raw Data
Handle
20406038
Attached To
rLAMMPS lammps
fix_rx.cpp
View Options
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fix_rx.h"
#include "atom.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "force.h"
#include "memory.h"
#include "comm.h"
#include "update.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_special.h"
#include "pair_dpd_fdt_energy.h"
#include <float.h> // DBL_EPSILON
#include <vector> // std::vector<>
#include <algorithm> // std::max
#include <cmath> // std::fmod
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathSpecial;
enum{NONE,HARMONIC};
enum{LUCY};
#define MAXLINE 1024
#define DELTA 4
#define SparseKinetics_enableIntegralReactions (true)
#define SparseKinetics_invalidIndex (-1)
namespace /* anonymous */
{
typedef double TimerType;
TimerType getTimeStamp(void) { return MPI_Wtime(); }
double getElapsedTime( const TimerType &t0, const TimerType &t1) { return t1-t0; }
} // end namespace
/* ---------------------------------------------------------------------- */
FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7 || narg > 12) error->all(FLERR,"Illegal fix rx command");
restart_peratom = 1;
nevery = 1;
nreactions = maxparam = 0;
params = NULL;
mol2param = NULL;
pairDPDE = NULL;
id_fix_species = NULL;
id_fix_species_old = NULL;
const int Verbosity = 1;
// Keep track of the argument list.
int iarg = 3;
// Read the kinetic file in arg[3].
kineticsFile = arg[iarg++];
// Determine the local temperature averaging method in arg[4].
wtFlag = 0;
localTempFlag = NONE;
{
char *word = arg[iarg++];
if (strcmp(word,"none") == 0){
wtFlag = 0;
localTempFlag = NONE;
}
else if (strcmp(word,"lucy") == 0){
wtFlag = LUCY;
localTempFlag = HARMONIC;
}
else
error->all(FLERR,"Illegal fix rx local temperature weighting technique");
}
// Select either sparse and dense matrix
// representations of the stoichiometric matrix.
useSparseKinetics = true;
{
char *word = arg[iarg++];
if (strcmp(word,"sparse") == 0)
useSparseKinetics = true;
else if (strcmp(word,"dense") == 0)
useSparseKinetics = false;
else {
std::string errmsg = "Illegal command " + std::string(word)
+ " expected \"sparse\" or \"dense\"\n";
error->all(FLERR, errmsg.c_str());
}
if (comm->me == 0 and Verbosity > 1){
std::string msg = "FixRX: matrix format is ";
if (useSparseKinetics)
msg += std::string("sparse");
else
msg += std::string("dense");
error->message(FLERR, msg.c_str());
}
}
// Determine the ODE solver/stepper strategy in arg[6].
odeIntegrationFlag = ODE_LAMMPS_RK4;
{
char *word = arg[iarg++];
if (strcmp(word,"lammps_rk4") == 0 || strcmp(word,"rk4") == 0)
odeIntegrationFlag = ODE_LAMMPS_RK4;
else if (strcmp(word,"lammps_rkf45") == 0 || strcmp(word,"rkf45") == 0)
odeIntegrationFlag = ODE_LAMMPS_RKF45;
else {
std::string errmsg = "Illegal ODE integration type: " + std::string(word);
error->all(FLERR, errmsg.c_str());
}
}
/// Set the default ODE parameters here. Modify with arg[].
/// 'minSteps' has a different meaning for RK4 and RKF45.
/// RK4: This is the # of steps that will be taken with h = dt_dpd / minSteps;
/// RKF45: This sets as h0 = dt_dpd / minSteps. If minSteps == 0, RKF45 will
/// estimate h0 internally. h will be adjusted as needed on subsequent steps.
minSteps = 1;
maxIters = 100;
relTol = 1.0e-6;
absTol = 1.0e-8;
diagnosticFrequency = 0;
for (int i = 0; i < numDiagnosticCounters; ++i){
diagnosticCounter[i] = 0;
diagnosticCounterPerODE[i] = NULL;
}
if (odeIntegrationFlag == ODE_LAMMPS_RK4 && narg==8){
char *word = arg[iarg++];
minSteps = atoi( word );
if (comm->me == 0 and Verbosity > 1){
char msg[128];
sprintf(msg, "FixRX: RK4 numSteps= %d", minSteps);
error->message(FLERR, msg);
}
}
else if (odeIntegrationFlag == ODE_LAMMPS_RK4 && narg>8){
error->all(FLERR,"Illegal fix rx command. Too many arguments for RK4 solver.");
}
else if (odeIntegrationFlag == ODE_LAMMPS_RKF45){
// Must have four options.
if (narg < 11)
error->all(FLERR,"Illegal fix rx command. Too few arguments for RKF45 solver.");
minSteps = atoi( arg[iarg++] );
maxIters = atoi( arg[iarg++] );
relTol = strtod( arg[iarg++], NULL);
absTol = strtod( arg[iarg++], NULL);
if (iarg < narg)
diagnosticFrequency = atoi( arg[iarg++] );
// maxIters must be at least minSteps.
maxIters = std::max( minSteps, maxIters );
if (comm->me == 0 and Verbosity > 1){
//printf("FixRX: RKF45 minSteps= %d maxIters= %d absTol= %e relTol= %e\n", minSteps, maxIters, absTol, relTol);
char msg[128];
sprintf(msg, "FixRX: RKF45 minSteps= %d maxIters= %d relTol= %.1e absTol= %.1e diagnosticFrequency= %d", minSteps, maxIters, relTol, absTol, diagnosticFrequency);
error->message(FLERR, msg);
}
}
// Initialize/Create the sparse matrix database.
sparseKinetics_nu = NULL;
sparseKinetics_nuk = NULL;
sparseKinetics_inu = NULL;
sparseKinetics_isIntegralReaction = NULL;
sparseKinetics_maxReactants = 0;
sparseKinetics_maxProducts = 0;
sparseKinetics_maxSpecies = 0;
}
/* ---------------------------------------------------------------------- */
FixRX::~FixRX()
{
// De-Allocate memory to prevent memory leak
for (int ii = 0; ii < nreactions; ii++){
delete [] stoich[ii];
delete [] stoichReactants[ii];
delete [] stoichProducts[ii];
}
delete [] Arr;
delete [] nArr;
delete [] Ea;
delete [] tempExp;
delete [] stoich;
delete [] stoichReactants;
delete [] stoichProducts;
delete [] kR;
delete [] id_fix_species;
delete [] id_fix_species_old;
if (useSparseKinetics){
memory->destroy( sparseKinetics_nu );
memory->destroy( sparseKinetics_nuk );
memory->destroy( sparseKinetics_inu );
memory->destroy( sparseKinetics_isIntegralReaction );
}
}
/* ---------------------------------------------------------------------- */
void FixRX::post_constructor()
{
int maxspecies = 1000;
int nUniqueSpecies = 0;
bool match;
for (int i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
error->all(FLERR,"fix rx cannot be combined with fix property/atom");
char **tmpspecies = new char*[maxspecies];
for(int jj=0; jj < maxspecies; jj++)
tmpspecies[jj] = NULL;
// open file on proc 0
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = force->open_potential(kineticsFile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",kineticsFile);
error->one(FLERR,str);
}
}
// Assign species names to tmpspecies array and determine the number of unique species
int n,nwords;
char line[MAXLINE],*ptr;
int eof = 0;
char * word;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// words = ptrs to all words in line
nwords = 0;
word = strtok(line," \t\n\r\f");
while (word != NULL){
word = strtok(NULL, " \t\n\r\f");
match=false;
for(int jj=0;jj<nUniqueSpecies;jj++){
if(strcmp(word,tmpspecies[jj])==0){
match=true;
break;
}
}
if(!match){
if(nUniqueSpecies+1>=maxspecies)
error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
strcpy(tmpspecies[nUniqueSpecies],word);
nUniqueSpecies++;
}
word = strtok(NULL, " \t\n\r\f");
if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) break;
word = strtok(NULL, " \t\n\r\f");
}
}
atom->nspecies_dpd = nUniqueSpecies;
nspecies = atom->nspecies_dpd;
// new id = fix-ID + FIX_STORE_ATTRIBUTE
// new fix group = group for this fix
id_fix_species = NULL;
id_fix_species_old = NULL;
n = strlen(id) + strlen("_SPECIES") + 1;
id_fix_species = new char[n];
n = strlen(id) + strlen("_SPECIES_OLD") + 1;
id_fix_species_old = new char[n];
strcpy(id_fix_species,id);
strcat(id_fix_species,"_SPECIES");
strcpy(id_fix_species_old,id);
strcat(id_fix_species_old,"_SPECIES_OLD");
char **newarg = new char*[nspecies+5];
char **newarg2 = new char*[nspecies+5];
newarg[0] = id_fix_species;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "property/atom";
newarg2[0] = id_fix_species_old;
newarg2[1] = group->names[igroup];
newarg2[2] = (char *) "property/atom";
for(int ii=0; ii<nspecies; ii++){
char str1[2+strlen(tmpspecies[ii])+1];
char str2[2+strlen(tmpspecies[ii])+4];
strcpy(str1,"d_");
strcpy(str2,"d_");
strncat(str1,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,"Old",3);
newarg[ii+3] = new char[strlen(str1)+1];
newarg2[ii+3] = new char[strlen(str2)+1];
strcpy(newarg[ii+3],str1);
strcpy(newarg2[ii+3],str2);
}
newarg[nspecies+3] = (char *) "ghost";
newarg[nspecies+4] = (char *) "yes";
newarg2[nspecies+3] = (char *) "ghost";
newarg2[nspecies+4] = (char *) "yes";
modify->add_fix(nspecies+5,newarg);
fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
modify->add_fix(nspecies+5,newarg2);
fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];
if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
for(int jj=0;jj<nspecies;jj++) {
delete[] tmpspecies[jj];
delete[] newarg[jj+3];
delete[] newarg2[jj+3];
}
delete[] newarg;
delete[] newarg2;
delete[] tmpspecies;
read_file( kineticsFile );
if (useSparseKinetics)
this->initSparse();
// set comm size needed by this Pair
comm_forward = nspecies*2;
comm_reverse = 2;
}
/* ---------------------------------------------------------------------- */
void FixRX::initSparse()
{
const int Verbosity = 1;
if (comm->me == 0 and Verbosity > 1){
for (int k = 0; k < nspecies; ++k)
printf("atom->dname[%d]= %s\n", k, atom->dname[k]);
printf("stoich[][]\n");
for (int i = 0; i < nreactions; ++i){
int nreac_i = 0, nprod_i = 0;
printf("%d: ", i);
for (int k = 0; k < nspecies; ++k){
printf(" %g", stoich[i][k]);
if (stoich[i][k] < 0.0) nreac_i++;
else if (stoich[i][k] > 0.0) nprod_i++;
}
printf(" : %d %d\n", nreac_i, nprod_i);
}
printf("stoichReactants[][]\n");
for (int i = 0; i < nreactions; ++i){
int nreac_i = 0;
printf("%d: ", i);
for (int k = 0; k < nspecies; ++k){
printf(" %g", stoichReactants[i][k]);
if (stoichReactants[i][k] > 0.0) nreac_i++;
}
printf(" : %d\n", nreac_i);
}
printf("stoichProducts[][]\n");
for (int i = 0; i < nreactions; ++i){
int nprod_i = 0;
printf("%d: ", i);
for (int k = 0; k < nspecies; ++k){
printf(" %g", stoichProducts[i][k]);
if (stoichProducts[i][k] > 0.0) nprod_i++;
}
printf(" : %d\n", nprod_i);
}
} // if (Verbose)
// 1) Measure the sparsity of stoich[][]
int nzeros = 0;
int mxprod = 0;
int mxreac = 0;
int mxspec = 0;
int nIntegral = 0;
for (int i = 0; i < nreactions; ++i){
int nreac_i = 0, nprod_i = 0;
std::string pstr, rstr;
bool allAreIntegral = true;
for (int k = 0; k < nspecies; ++k){
if (stoichReactants[i][k] == 0 and stoichProducts[i][k] == 0)
nzeros++;
if (stoichReactants[i][k] > 0.0){
allAreIntegral &= (std::fmod( stoichReactants[i][k], 1.0 ) == 0.0);
nreac_i++;
if (rstr.length() > 0)
rstr += " + ";
char digit[6];
sprintf(digit, "%4.1f ", stoichReactants[i][k]); rstr += digit;
rstr += atom->dname[k];
}
if (stoichProducts[i][k] > 0.0){
allAreIntegral &= (std::fmod( stoichProducts[i][k], 1.0 ) == 0.0);
nprod_i++;
if (pstr.length() > 0)
pstr += " + ";
char digit[6];
sprintf(digit, "%4.1f ", stoichProducts[i][k]); pstr += digit;
pstr += atom->dname[k];
}
}
if (comm->me == 0 and Verbosity > 1)
printf("rx%3d: %d %d %d ... %s %s %s\n", i, nreac_i, nprod_i, allAreIntegral, rstr.c_str(), /*reversible[i]*/ (false) ? "<=>" : "=", pstr.c_str());
mxreac = std::max( mxreac, nreac_i );
mxprod = std::max( mxprod, nprod_i );
mxspec = std::max( mxspec, nreac_i + nprod_i );
if (allAreIntegral) nIntegral++;
}
if (comm->me == 0 and Verbosity > 1){
char msg[256];
sprintf(msg, "FixRX: Sparsity of Stoichiometric Matrix= %.1f%% non-zeros= %d nspecies= %d nreactions= %d maxReactants= %d maxProducts= %d maxSpecies= %d integralReactions= %d", 100*(double(nzeros) / (nspecies * nreactions)), nzeros, nspecies, nreactions, mxreac, mxprod, (mxreac + mxprod), SparseKinetics_enableIntegralReactions);
error->message(FLERR, msg);
}
// Allocate the sparse matrix data.
{
sparseKinetics_maxSpecies = (mxreac + mxprod);
sparseKinetics_maxReactants = mxreac;
sparseKinetics_maxProducts = mxprod;
memory->create( sparseKinetics_nu , nreactions, sparseKinetics_maxSpecies, "sparseKinetics_nu");
memory->create( sparseKinetics_nuk, nreactions, sparseKinetics_maxSpecies, "sparseKinetics_nuk");
for (int i = 0; i < nreactions; ++i)
for (int k = 0; k < sparseKinetics_maxSpecies; ++k){
sparseKinetics_nu [i][k] = 0.0;
sparseKinetics_nuk[i][k] = SparseKinetics_invalidIndex; // Initialize with an invalid index.
}
if (SparseKinetics_enableIntegralReactions){
memory->create( sparseKinetics_inu, nreactions, sparseKinetics_maxSpecies, "sparseKinetics_inu");
memory->create( sparseKinetics_isIntegralReaction, nreactions, "sparseKinetics_isIntegralReaction");
for (int i = 0; i < nreactions; ++i){
sparseKinetics_isIntegralReaction[i] = false;
for (int k = 0; k < sparseKinetics_maxSpecies; ++k)
sparseKinetics_inu[i][k] = 0;
}
}
}
// Measure the distribution of the # of moles for the ::fastpowi function.
std::vector<int> nu_bin(10);
for (int i = 0; i < nreactions; ++i){
int nreac_i = 0, nprod_i = 0;
bool isIntegral_i = true;
for (int k = 0; k < nspecies; ++k){
if (stoichReactants[i][k] > 0.0){
const int idx = nreac_i;
sparseKinetics_nu [i][idx] = stoichReactants[i][k];
sparseKinetics_nuk[i][idx] = k;
isIntegral_i &= (std::fmod( stoichReactants[i][k], 1.0 ) == 0.0);
if (SparseKinetics_enableIntegralReactions){
sparseKinetics_inu[i][idx] = (int)sparseKinetics_nu[i][idx];
if (isIntegral_i){
if (sparseKinetics_inu[i][idx] >= nu_bin.size())
nu_bin.resize( sparseKinetics_inu[i][idx] );
nu_bin[ sparseKinetics_inu[i][idx] ] ++;
}
}
nreac_i++;
}
if (stoichProducts[i][k] > 0.0){
const int idx = sparseKinetics_maxReactants + nprod_i;
sparseKinetics_nu [i][idx] = stoichProducts[i][k];
sparseKinetics_nuk[i][idx] = k;
isIntegral_i &= (std::fmod( sparseKinetics_nu[i][idx], 1.0 ) == 0.0);
if (SparseKinetics_enableIntegralReactions){
sparseKinetics_inu[i][idx] = (int) sparseKinetics_nu[i][idx];
if (isIntegral_i){
if (sparseKinetics_inu[i][idx] >= nu_bin.size())
nu_bin.resize( sparseKinetics_inu[i][idx] );
nu_bin[ sparseKinetics_inu[i][idx] ] ++;
}
}
nprod_i++;
}
}
if (SparseKinetics_enableIntegralReactions)
sparseKinetics_isIntegralReaction[i] = isIntegral_i;
}
if (comm->me == 0 and Verbosity > 1){
for (int i = 1; i < nu_bin.size(); ++i)
if (nu_bin[i] > 0)
printf("nu_bin[%d] = %d\n", i, nu_bin[i]);
for (int i = 0; i < nreactions; ++i){
std::string pstr, rstr;
for (int kk = 0; kk < sparseKinetics_maxReactants; kk++){
const int k = sparseKinetics_nuk[i][kk];
if (k != SparseKinetics_invalidIndex){
if (rstr.length() > 0)
rstr += " + ";
char digit[6];
if (SparseKinetics_enableIntegralReactions and sparseKinetics_isIntegralReaction[i])
sprintf(digit,"%d ", sparseKinetics_inu[i][kk]);
else
sprintf(digit,"%4.1f ", sparseKinetics_nu[i][kk]);
rstr += digit;
rstr += atom->dname[k];
}
}
for (int kk = sparseKinetics_maxReactants; kk < sparseKinetics_maxSpecies; kk++){
const int k = sparseKinetics_nuk[i][kk];
if (k != SparseKinetics_invalidIndex){
if (pstr.length() > 0)
pstr += " + ";
char digit[6];
if (SparseKinetics_enableIntegralReactions and sparseKinetics_isIntegralReaction[i])
sprintf(digit,"%d ", sparseKinetics_inu[i][kk]);
else
sprintf(digit,"%4.1f ", sparseKinetics_nu[i][kk]);
pstr += digit;
pstr += atom->dname[k];
}
}
if (comm->me == 0 and Verbosity > 1)
printf("rx%3d: %s %s %s\n", i, rstr.c_str(), /*reversible[i]*/ (false) ? "<=>" : "=", pstr.c_str());
}
// end for nreactions
}
// end if Verbose
}
/* ---------------------------------------------------------------------- */
int FixRX::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRX::init()
{
pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1);
if (pairDPDE == NULL)
error->all(FLERR,"Must use pair_style dpd/fdt/energy with fix rx");
bool eos_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"eos/table/rx") == 0) eos_flag = true;
if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");
}
/* ---------------------------------------------------------------------- */
void FixRX::setup_pre_force(int vflag)
{
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int *mask = atom->mask;
int newton_pair = force->newton_pair;
double tmp;
int ii;
if(localTempFlag){
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
}
for (int id = 0; id < nlocal; id++)
for (int ispecies=0; ispecies<nspecies; ispecies++){
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;
if (odeIntegrationFlag == ODE_LAMMPS_RK4)
rk4(i,NULL);
else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
rkf45(i,NULL);
}
// Communicate the updated momenta and velocities to all nodes
comm->forward_comm_fix(this);
if(localTempFlag) delete [] dpdThetaLocal;
}
/* ---------------------------------------------------------------------- */
void FixRX::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;
int ii;
double theta;
if(localTempFlag){
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
}
TimerType timer_localTemperature = getTimeStamp();
// Zero the counters for the ODE solvers.
this->nSteps = this->nIters = this->nFuncs = this->nFails = 0;
if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency == 1)
{
memory->create( diagnosticCounterPerODE[StepSum], nlocal, "FixRX::diagnosticCounterPerODE");
memory->create( diagnosticCounterPerODE[FuncSum], nlocal, "FixRX::diagnosticCounterPerODE");
}
double *rwork = new double[8*nspecies + nreactions];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
if (localTempFlag)
theta = dpdThetaLocal[i];
else
theta = dpdTheta[i];
//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);
if (odeIntegrationFlag == ODE_LAMMPS_RK4)
rk4(i,rwork);
else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
rkf45(i,rwork);
}
TimerType timer_ODE = getTimeStamp();
delete [] rwork;
// Communicate the updated momenta and velocities to all nodes
comm->forward_comm_fix(this);
if(localTempFlag) delete [] dpdThetaLocal;
double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE);
// Warn the user if a failure was detected in the ODE solver.
if (nFails > 0){
char sbuf[128];
sprintf(sbuf,"in FixRX::pre_force, ODE solver failed for %d atoms.", nFails);
error->warning(FLERR, sbuf);
}
// Compute and report ODE diagnostics, if requested.
if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency != 0){
// Update the counters.
diagnosticCounter[StepSum] += nSteps;
diagnosticCounter[FuncSum] += nFuncs;
diagnosticCounter[TimeSum] += time_ODE;
diagnosticCounter[AtomSum] += nlocal;
diagnosticCounter[numDiagnosticCounters-1] ++;
if ( (diagnosticFrequency > 0 &&
((update->ntimestep - update->firststep) % diagnosticFrequency) == 0) ||
(diagnosticFrequency < 0 && update->ntimestep == update->laststep) )
this->odeDiagnostics();
for (int i = 0; i < numDiagnosticCounters; ++i)
if (diagnosticCounterPerODE[i])
memory->destroy( diagnosticCounterPerODE[i] );
}
}
/* ---------------------------------------------------------------------- */
void FixRX::read_file(char *file)
{
nreactions = 0;
// open file on proc 0
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",file);
error->one(FLERR,str);
}
}
// Count the number of reactions from kinetics file
int n,nwords,ispecies;
char line[MAXLINE],*ptr;
int eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
nreactions++;
}
// open file on proc 0
if (comm->me == 0) fp = force->open_potential(file);
// read each reaction from kinetics file
eof=0;
char * word;
double tmpStoich;
double sign;
Arr = new double[nreactions];
nArr = new double[nreactions];
Ea = new double[nreactions];
tempExp = new double[nreactions];
stoich = new double*[nreactions];
stoichReactants = new double*[nreactions];
stoichProducts = new double*[nreactions];
for (int ii=0;ii<nreactions;ii++){
stoich[ii] = new double[nspecies];
stoichReactants[ii] = new double[nspecies];
stoichProducts[ii] = new double[nspecies];
}
kR = new double[nreactions];
for (int ii=0;ii<nreactions;ii++){
for (int jj=0;jj<nspecies;jj++){
stoich[ii][jj] = 0.0;
stoichReactants[ii][jj] = 0.0;
stoichProducts[ii][jj] = 0.0;
}
}
nreactions=0;
sign = -1.0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// words = ptrs to all words in line
nwords = 0;
word = strtok(line," \t\n\r\f");
while (word != NULL){
tmpStoich = atof(word);
word = strtok(NULL, " \t\n\r\f");
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(word,&atom->dname[ispecies][0]) == 0){
stoich[nreactions][ispecies] += sign*tmpStoich;
if(sign<0.0)
stoichReactants[nreactions][ispecies] += tmpStoich;
else stoichProducts[nreactions][ispecies] += tmpStoich;
break;
}
}
if(ispecies==nspecies){
if (comm->me) {
fprintf(stderr,"%s mol fraction is not found in data file\n",word);
fprintf(stderr,"nspecies=%d ispecies=%d\n",nspecies,ispecies);
}
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");
Arr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
nArr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
Ea[nreactions] = atof(word);
sign = -1.0;
break;
}
word = strtok(NULL, " \t\n\r\f");
}
nreactions++;
}
}
/* ---------------------------------------------------------------------- */
void FixRX::setupParams()
{
int i,j,n;
// set mol2param for all combinations
// must be a single exact match to lines read from file
memory->destroy(mol2param);
memory->create(mol2param,nspecies,"pair:mol2param");
for (i = 0; i < nspecies; i++) {
n = -1;
for (j = 0; j < nreactions; j++) {
if (i == params[j].ispecies) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = j;
}
}
mol2param[i] = n;
}
}
/* ---------------------------------------------------------------------- */
void FixRX::rk4(int id, double *rwork)
{
double *k1 = NULL;
if (rwork == NULL)
k1 = new double[6*nspecies + nreactions];
else
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);
// Update ConcOld
for (int ispecies = 0; ispecies < nspecies; ispecies++)
{
const double tmp = atom->dvector[ispecies][id];
atom->dvector[ispecies+nspecies][id] = tmp;
y[ispecies] = tmp;
}
// Run the requested steps with h.
for (int step = 0; step < numSteps; step++)
{
// k1
rhs(0.0,y,k1,dummyArray);
// k2
for (int ispecies = 0; ispecies < nspecies; ispecies++)
yp[ispecies] = y[ispecies] + 0.5*h*k1[ispecies];
rhs(0.0,yp,k2,dummyArray);
// k3
for (int ispecies = 0; ispecies < nspecies; ispecies++)
yp[ispecies] = y[ispecies] + 0.5*h*k2[ispecies];
rhs(0.0,yp,k3,dummyArray);
// k4
for (int ispecies = 0; ispecies < nspecies; ispecies++)
yp[ispecies] = y[ispecies] + h*k3[ispecies];
rhs(0.0,yp,k4,dummyArray);
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);
} // end for (int step...
// Store the solution back in atom->dvector.
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if(y[ispecies] < -1.0e-10)
error->one(FLERR,"Computed concentration in RK4 solver is < -1.0e-10");
else if(y[ispecies] < 1e-15)
y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
if (rwork == NULL)
delete [] k1;
}
/* ---------------------------------------------------------------------- */
// f1 = dt*f(t,x)
// f2 = dt*f(t+ c20*dt,x + c21*f1)
// f3 = dt*f(t+ c30*dt,x + c31*f1 + c32*f2)
// f4 = dt*f(t+ c40*dt,x + c41*f1 + c42*f2 + c43*f3)
// f5 = dt*f(t+dt,x + c51*f1 + c52*f2 + c53*f3 + c54*f4)
// f6 = dt*f(t+ c60*dt,x + c61*f1 + c62*f2 + c63*f3 + c64*f4 + c65*f5)
//
// fifth-order runge-kutta integration
// x5 = x + b1*f1 + b3*f3 + b4*f4 + b5*f5 + b6*f6
// fourth-order runge-kutta integration
// x = x + a1*f1 + a3*f3 + a4*f4 + a5*f5
void FixRX::rkf45_step (const int neq, const double h, double y[], double y_out[], double rwk[], void* v_param)
{
const double c21=0.25;
const double c31=0.09375;
const double c32=0.28125;
const double c41=0.87938097405553;
const double c42=-3.2771961766045;
const double c43=3.3208921256258;
const double c51=2.0324074074074;
const double c52=-8.0;
const double c53=7.1734892787524;
const double c54=-0.20589668615984;
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 a3=0.54892787524366;
const double a4=0.5353313840156;
const double a5=-0.2;
const double b1=0.11851851851852;
const double b3=0.51898635477583;
const double b4=0.50613149034201;
const double b5=-0.18;
const double b6=0.036363636363636;
// local dependent variables (5 total)
double* f1 = &rwk[ 0];
double* f2 = &rwk[ neq];
double* f3 = &rwk[2*neq];
double* f4 = &rwk[3*neq];
double* f5 = &rwk[4*neq];
double* f6 = &rwk[5*neq];
// scratch for the intermediate solution.
//double* ytmp = &rwk[6*neq];
double* ytmp = y_out;
// 1)
rhs (0.0, y, f1, v_param);
for (int k = 0; k < neq; k++){
f1[k] *= h;
ytmp[k] = y[k] + c21 * f1[k];
}
// 2)
rhs(0.0, ytmp, f2, v_param);
for (int k = 0; k < neq; k++){
f2[k] *= h;
ytmp[k] = y[k] + c31 * f1[k] + c32 * f2[k];
}
// 3)
rhs(0.0, ytmp, f3, v_param);
for (int k = 0; k < neq; k++) {
f3[k] *= h;
ytmp[k] = y[k] + c41 * f1[k] + c42 * f2[k] + c43 * f3[k];
}
// 4)
rhs(0.0, ytmp, f4, v_param);
for (int k = 0; k < neq; k++) {
f4[k] *= h;
ytmp[k] = y[k] + c51 * f1[k] + c52 * f2[k] + c53 * f3[k] + c54 * f4[k];
}
// 5)
rhs(0.0, ytmp, f5, v_param);
for (int k = 0; k < neq; k++) {
f5[k] *= h;
ytmp[k] = y[k] + c61*f1[k] + c62*f2[k] + c63*f3[k] + c64*f4[k] + c65*f5[k];
}
// 6)
rhs(0.0, ytmp, f6, v_param);
for (int k = 0; k < neq; k++)
{
//const double f6 = h * ydot[k];
f6[k] *= h;
// 5th-order solution.
const double r5 = b1*f1[k] + b3*f3[k] + b4*f4[k] + b5*f5[k] + b6*f6[k];
// 4th-order solution.
const double r4 = a1*f1[k] + a3*f3[k] + a4*f4[k] + a5*f5[k];
// Truncation error: difference between 4th and 5th-order solutions.
rwk[k] = fabs(r5 - r4);
// Update solution.
//y_out[k] = y[k] + r5; // Local extrapolation
y_out[k] = y[k] + r4;
}
return;
}
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)
{
// 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.
// Adjust upper bound based on ydot ...
double hg = sqrt(hmin*hmax);
//if (hmax < hmin)
//{
// h0 = hg;
// return;
//}
// Start iteration to find solution to ... {WRMS norm of (h0^2 y'' / 2)} = 1
double *ydot = rwk;
double *y1 = ydot + neq;
double *ydot1 = y1 + neq;
const int max_iters = 10;
bool hnew_is_ok = false;
double hnew = hg;
int iter = 0;
// compute ydot at t=t0
rhs (t, y, ydot, v_params);
while(1)
{
// Estimate y'' with finite-difference ...
for (int k = 0; k < neq; k++)
y1[k] = y[k] + hg * ydot[k];
// compute y' at t1
rhs (t + hg, y1, ydot1, v_params);
// Compute WRMS norm of y''
double yddnrm = 0.0;
for (int k = 0; k < neq; k++){
double ydd = (ydot1[k] - ydot[k]) / hg;
double wterr = ydd / (relTol * fabs( y[k] ) + absTol);
yddnrm += wterr * wterr;
}
yddnrm = sqrt( yddnrm / double(neq) );
//std::cout << "iter " << _iter << " hg " << hg << " y'' " << yddnrm << std::endl;
//std::cout << "ydot " << ydot[neq-1] << std::endl;
// should we accept this?
if (hnew_is_ok || iter == max_iters){
hnew = hg;
if (iter == max_iters)
fprintf(stderr, "ERROR_HIN_MAX_ITERS\n");
break;
}
// Get the new value of h ...
hnew = (yddnrm*hmax*hmax > 2.0) ? sqrt(2.0 / yddnrm) : sqrt(hg * hmax);
// test the stopping conditions.
double hrat = hnew / hg;
// Accept this value ... the bias factor should bring it within range.
if ( (hrat > 0.5) && (hrat < 2.0) )
hnew_is_ok = true;
// If y'' is still bad after a few iterations, just accept h and give up.
if ( (iter > 1) && hrat > 2.0 ) {
hnew = hg;
hnew_is_ok = true;
}
//printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax);
hg = hnew;
iter ++;
}
// bound and bias estimate
h0 = hnew * 0.5;
h0 = fmax(h0, hmin);
h0 = fmin(h0, hmax);
//printf("h0=%e, hmin=%e, hmax=%e\n", h0, hmin, hmax);
return (iter + 1);
}
void FixRX::odeDiagnostics(void)
{
TimerType timer_start = getTimeStamp();
// Compute:
// 1) Average # of ODE integrator steps and RHS evaluations per atom globally.
// 2) RMS # of ...
// 3) Average # of ODE steps and RHS evaluations per MPI task.
// 4) RMS # of ODE steps and RHS evaluations per MPI task.
// 5) MAX # of ODE steps and RHS evaluations per MPI task.
//
// ... 1,2 are for ODE control diagnostics.
// ... 3-5 are for load balancing diagnostics.
//
// To do this, we'll need to
// a) Allreduce (sum) the sum of nSteps / nFuncs. Dividing by atom->natoms
// gives the avg # of steps/funcs per atom globally.
// b) Reduce (sum) to root the sum of squares of the differences.
// i) Sum_i (steps_i - avg_steps_global)^2
// ii) Sum_i (funcs_i - avg_funcs_global)^2
// iii) (avg_steps_local - avg_steps_global)^2
// iv) (avg_funcs_local - avg_funcs_global)^2
const int numCounters = numDiagnosticCounters-1;
// # of time-steps for averaging.
const int nTimes = this->diagnosticCounter[numDiagnosticCounters-1];
// # of ODE's per time-step (on average).
//const int nODEs = this->diagnosticCounter[AtomSum] / nTimes;
// Sum up the sums from each task.
double sums[numCounters];
double my_vals[numCounters];
double max_per_proc[numCounters];
double min_per_proc[numCounters];
// Compute counters per dpd time-step.
for (int i = 0; i < numCounters; ++i){
my_vals[i] = this->diagnosticCounter[i] / nTimes;
//printf("my sum[%d] = %f %d\n", i, my_vals[i], comm->me);
}
MPI_Allreduce (my_vals, sums, numCounters, MPI_DOUBLE, MPI_SUM, world);
MPI_Reduce (my_vals, max_per_proc, numCounters, MPI_DOUBLE, MPI_MAX, 0, world);
MPI_Reduce (my_vals, min_per_proc, numCounters, MPI_DOUBLE, MPI_MIN, 0, world);
const double nODEs = sums[numCounters-1];
double avg_per_atom[numCounters], avg_per_proc[numCounters];
// Averages per-ODE and per-proc per time-step.
for (int i = 0; i < numCounters; ++i){
avg_per_atom[i] = sums[i] / nODEs;
avg_per_proc[i] = sums[i] / comm->nprocs;
}
// Sum up the differences from each task.
double sum_sq[2*numCounters];
double my_sum_sq[2*numCounters];
for (int i = 0; i < numCounters; ++i){
double diff_i = my_vals[i] - avg_per_proc[i];
my_sum_sq[i] = diff_i * diff_i;
}
double max_per_ODE[numCounters], min_per_ODE[numCounters];
// Process the per-ODE RMS of the # of steps/funcs
if (diagnosticFrequency == 1){
double my_max[numCounters], my_min[numCounters];
const int nlocal = atom->nlocal;
const int *mask = atom->mask;
for (int i = 0; i < numCounters; ++i){
my_sum_sq[i+numCounters] = 0;
my_max[i] = 0;
my_min[i] = DBL_MAX;
if (diagnosticCounterPerODE[i] != NULL){
for (int j = 0; j < nlocal; ++j)
if (mask[j] & groupbit){
double diff = double(diagnosticCounterPerODE[i][j]) - avg_per_atom[i];
my_sum_sq[i+numCounters] += diff*diff;
my_max[i] = std::max( my_max[i], (double)diagnosticCounterPerODE[i][j] );
my_min[i] = std::min( my_min[i], (double)diagnosticCounterPerODE[i][j] );
}
}
}
MPI_Reduce (my_sum_sq, sum_sq, 2*numCounters, MPI_DOUBLE, MPI_SUM, 0, world);
MPI_Reduce (my_max, max_per_ODE, numCounters, MPI_DOUBLE, MPI_MAX, 0, world);
MPI_Reduce (my_min, min_per_ODE, numCounters, MPI_DOUBLE, MPI_MIN, 0, world);
}
else
MPI_Reduce (my_sum_sq, sum_sq, numCounters, MPI_DOUBLE, MPI_SUM, 0, world);
TimerType timer_stop = getTimeStamp();
double time_local = getElapsedTime( timer_start, timer_stop );
if (comm->me == 0){
char smesg[128];
#define print_mesg(smesg) {\
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)");
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]);
print_mesg(smesg);
// only valid for single time-step!
if (diagnosticFrequency == 1){
double rms_per_ODE[numCounters];
for (int i = 0; i < numCounters; ++i)
rms_per_ODE[i] = sqrt( sum_sq[i+numCounters] / nODEs );
sprintf(smesg, " RMS per ODE : %-12.5g | %-12.5g ", rms_per_ODE[0], rms_per_ODE[1]);
print_mesg(smesg);
sprintf(smesg, " MAX per ODE : %-12.5g | %-12.5g ", max_per_ODE[0], max_per_ODE[1]);
print_mesg(smesg);
sprintf(smesg, " MIN per ODE : %-12.5g | %-12.5g ", min_per_ODE[0], min_per_ODE[1]);
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]);
print_mesg(smesg);
if (comm->nprocs > 1){
double rms_per_proc[numCounters];
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]);
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]);
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]);
print_mesg(smesg);
}
sprintf(smesg, " AVG'd over %d time-steps", nTimes);
print_mesg(smesg);
sprintf(smesg, " AVG'ing took %g sec", time_local);
print_mesg(smesg);
#undef print_mesg
}
// Reset the counters.
for (int i = 0; i < numDiagnosticCounters; ++i)
diagnosticCounter[i] = 0;
return;
}
void FixRX::rkf45(int id, double *rwork)
{
// Rounding coefficient.
const double uround = DBL_EPSILON;
// Adaption limit (shrink or grow)
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;
const int neq = nspecies;
// Update ConcOld and initialize the ODE solution vector y[].
for (int ispecies = 0; ispecies < nspecies; ispecies++){
const double tmp = atom->dvector[ispecies][id];
atom->dvector[ispecies+nspecies][id] = tmp;
y[ispecies] = tmp;
}
// Integration length.
const double t_stop = update->dt; // DPD time-step.
// Safety factor on the adaption. very specific but not necessary .. 0.9 is common.
const double hsafe = 0.840896415;
// Time rounding factor.
const double tround = t_stop * uround;
// Counters for diagnostics.
int nst = 0; // # of steps (accepted)
int nit = 0; // # of iterations total
int nfe = 0; // # of RHS evaluations
// Min/Max step-size limits.
const double h_min = 100.0 * tround;
const double h_max = (minSteps > 0) ? t_stop / double(minSteps) : t_stop;
// Set the initial step-size. 0 forces an internal estimate ... stable Euler step size.
double h = (minSteps > 0) ? t_stop / double(minSteps) : 0.0;
double t = 0.0;
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);
}
//printf("t= %e t_stop= %e h= %e\n", t, t_stop, h);
// Integrate until we reach the end time.
while (fabs(t - t_stop) > tround){
double *yout = y + neq;
double *eout = yout + neq;
// Take a trial step.
rkf45_step (neq, h, y, yout, eout, rhstmp);
// Estimate the solution error.
// ... weighted 2-norm of the error.
double err2 = 0.0;
for (int k = 0; k < neq; k++){
const double wterr = eout[k] / (relTol * fabs( y[k] ) + absTol);
err2 += wterr * wterr;
}
double err = fmax( uround, sqrt( err2 / double(nspecies) ));
// Accept the solution?
if (err <= 1.0 || h <= h_min){
t += h;
nst++;
for (int k = 0; k < neq; k++)
y[k] = yout[k];
}
// Adjust h for the next step.
double hfac = hsafe * sqrt( sqrt( 1.0 / err ) );
// Limit the adaption.
hfac = fmax( hfac, 1.0 / adaption_limit );
hfac = fmin( hfac, adaption_limit );
// Apply the adaption factor...
h *= hfac;
// Limit h.
h = fmin( h, h_max );
h = fmax( h, h_min );
// Stretch h if we're within 5% ... and we didn't just fail.
if (err <= 1.0 && (t + 1.05*h) > t_stop)
h = t_stop - t;
// And don't overshoot the end.
if (t + h > t_stop)
h = t_stop - t;
nit++;
nfe += 6;
if (maxIters && nit > maxIters){
//fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop);
nFails ++;
break;
// We should set an error here so that the solution is not used!
}
} // end while
nSteps += nst;
nIters += nit;
nFuncs += nfe;
//if (diagnosticFrequency == 1 && diagnosticCounterPerODE[StepSum] != NULL)
if (diagnosticCounterPerODE[StepSum] != NULL){
diagnosticCounterPerODE[StepSum][id] = nst;
diagnosticCounterPerODE[FuncSum][id] = nfe;
}
//printf("id= %d nst= %d nit= %d\n", id, nst, nit);
// Store the solution back in atom->dvector.
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if(y[ispecies] < -1.0e-10)
error->one(FLERR,"Computed concentration in RKF45 solver is < -1.0e-10");
else if(y[ispecies] < 1e-20)
y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
if (rwork == NULL)
delete [] y;
}
/* ---------------------------------------------------------------------- */
int FixRX::rhs(double t, const double *y, double *dydt, void *params)
{
// Use the sparse format instead.
if (useSparseKinetics)
return this->rhs_sparse( t, y, dydt, params);
else
return this->rhs_dense ( t, y, dydt, 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;
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];
for(int ispecies=0; ispecies<nspecies; ispecies++){
concentration = y[ispecies]/VDPD;
rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
}
rxnRateLaw[jrxn] = rxnRateLawForward;
}
// Construct the reaction rates for each species
for(int ispecies=0; ispecies<nspecies; ispecies++)
for(int jrxn=0; jrxn<nreactions; jrxn++)
dydt[ispecies] += stoich[jrxn][ispecies]*VDPD*rxnRateLaw[jrxn];
return 0;
}
/* ---------------------------------------------------------------------- */
int FixRX::rhs_sparse(double t, const double *y, double *dydt, void *v_params) const
{
double *_rxnRateLaw = (double *) v_params;
const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;
#define kFor (this->kR)
#define kRev (NULL)
#define rxnRateLaw (_rxnRateLaw)
#define conc (dydt)
#define maxReactants (this->sparseKinetics_maxReactants)
#define maxSpecies (this->sparseKinetics_maxSpecies)
#define nuk (this->sparseKinetics_nuk)
#define nu (this->sparseKinetics_nu)
#define inu (this->sparseKinetics_inu)
#define isIntegral(idx) (SparseKinetics_enableIntegralReactions \
&& this->sparseKinetics_isIntegralReaction[idx])
for (int k = 0; k < nspecies; ++k)
conc[k] = y[k] / VDPD;
// Construct the reaction rate laws
for (int i = 0; i < nreactions; ++i)
{
double rxnRateLawForward;
if (isIntegral(i)){
rxnRateLawForward = kFor[i] * powint( conc[ nuk[i][0] ], inu[i][0]);
for (int kk = 1; kk < maxReactants; ++kk){
const int k = nuk[i][kk];
if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= powint( conc[k], inu[i][kk] );
}
} else {
rxnRateLawForward = kFor[i] * pow( conc[ nuk[i][0] ], nu[i][0]);
for (int kk = 1; kk < maxReactants; ++kk){
const int k = nuk[i][kk];
if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= pow( conc[k], nu[i][kk] );
}
}
rxnRateLaw[i] = rxnRateLawForward;
}
// Construct the reaction rates for each species from the
// Stoichiometric matrix and ROP vector.
for (int k = 0; k < nspecies; ++k)
dydt[k] = 0.0;
for (int i = 0; i < nreactions; ++i){
// Reactants ...
dydt[ nuk[i][0] ] -= nu[i][0] * rxnRateLaw[i];
for (int kk = 1; kk < maxReactants; ++kk){
const int k = nuk[i][kk];
if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] -= nu[i][kk] * rxnRateLaw[i];
}
// Products ...
dydt[ nuk[i][maxReactants] ] += nu[i][maxReactants] * rxnRateLaw[i];
for (int kk = maxReactants+1; kk < maxSpecies; ++kk){
const int k = nuk[i][kk];
if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] += nu[i][kk] * rxnRateLaw[i];
}
}
// Add in the volume factor to convert to the proper units.
for (int k = 0; k < nspecies; ++k)
dydt[k] *= VDPD;
#undef kFor
#undef kRev
#undef rxnRateLaw
#undef conc
#undef maxReactants
#undef maxSpecies
#undef nuk
#undef nu
#undef inu
#undef isIntegral
//#undef invalidIndex
return 0;
}
/* ---------------------------------------------------------------------- */
void FixRX::computeLocalTemperature()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
int *type = atom->type;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int newton_pair = force->newton_pair;
// local temperature variables
double wij=0.0;
double *dpdTheta = atom->dpdTheta;
// Initialize the local density and local temperature arrays
if (newton_pair) {
sumWeights = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
sumWeights[ii] = 0.0;
} else {
sumWeights = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = 0.0;
}
inum = pairDPDE->list->inum;
ilist = pairDPDE->list->ilist;
numneigh = pairDPDE->list->numneigh;
firstneigh = pairDPDE->list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < pairDPDE->cutsq[itype][jtype]) {
double rcut = sqrt(pairDPDE->cutsq[itype][jtype]);
double rij = sqrt(rsq);
double ratio = rij/rcut;
// Lucy's Weight Function
if(wtFlag==LUCY){
wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio);
dpdThetaLocal[i] += wij/dpdTheta[j];
if (newton_pair || j < nlocal)
dpdThetaLocal[j] += wij/dpdTheta[i];
}
sumWeights[i] += wij;
if (newton_pair || j < nlocal)
sumWeights[j] += wij;
}
}
}
if (newton_pair) comm->reverse_comm_fix(this);
// self-interaction for local temperature
for (i = 0; i < nlocal; i++){
// Lucy Weight Function
if(wtFlag==LUCY){
wij = 1.0;
dpdThetaLocal[i] += wij / dpdTheta[i];
}
sumWeights[i] += wij;
// Normalized local temperature
dpdThetaLocal[i] = dpdThetaLocal[i] / sumWeights[i];
if(localTempFlag == HARMONIC)
dpdThetaLocal[i] = 1.0 / dpdThetaLocal[i];
}
delete [] sumWeights;
}
/* ---------------------------------------------------------------------- */
int FixRX::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int ii,jj,m;
double tmp;
m = 0;
for (ii = 0; ii < n; ii++) {
jj = list[ii];
for(int ispecies=0;ispecies<nspecies;ispecies++){
tmp = atom->dvector[ispecies][jj];
buf[m++] = tmp;
tmp = atom->dvector[ispecies+nspecies][jj];
buf[m++] = tmp;
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRX::unpack_forward_comm(int n, int first, double *buf)
{
int ii,m,last;
double tmp;
m = 0;
last = first + n ;
for (ii = first; ii < last; ii++){
for(int ispecies=0;ispecies<nspecies;ispecies++){
tmp = buf[m++];
atom->dvector[ispecies][ii] = tmp;
tmp = buf[m++];
atom->dvector[ispecies+nspecies][ii] = tmp;
}
}
}
/* ---------------------------------------------------------------------- */
int FixRX::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = dpdThetaLocal[i];
buf[m++] = sumWeights[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRX::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
dpdThetaLocal[j] += buf[m++];
sumWeights[j] += buf[m++];
}
}
Event Timeline
Log In to Comment