diff --git a/DropletCombustionTest7-Binary b/DropletCombustionTest7-Binary new file mode 100755 index 0000000..2ff96a9 Binary files /dev/null and b/DropletCombustionTest7-Binary differ diff --git a/Makefile b/Makefile index dd67f3d..0f5b8c9 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,7 @@ compiler =g++ CANTERA_DIR =/opt/scientific/cantera-2.4_gnu_blas IDA_DIR =/opt/scientific/sundials-3.1.1_intel_mkl -EXE =DropletCombustionTest7-binary +EXE =DropletCombustionTest7-C7_C12 #DESTDIR =~/bin DESTDIR = ../example diff --git a/UserData.cpp b/UserData.cpp index 85b5de6..8473e99 100644 --- a/UserData.cpp +++ b/UserData.cpp @@ -57,7 +57,14 @@ void freeUserData(UserData data){ fclose(data->timescaleOutput); printf("Characteristic Timescale Output File Cleared from Memory!\n"); } - + //if(data->rxnROPOutput!=NULL){ + // fclose(data->rxnROPOutput); + // printf("Reactions Rate of Progress Output File Cleared from Memory!\n"); + //} + //if(data->spROPOutput!=NULL){ + // fclose(data->spROPOutput); + // printf("Species Rate of Production Output File Cleared from Memory!\n"); + //} } free(data); /* Free the user data */ printf("\n\n"); @@ -259,7 +266,7 @@ UserData allocateUserData(FILE *input){ ier=parseNumber (input, "suppressAlg" , MAXBUFLEN, &data->suppressAlg); - if(data->setConstraints!=0 && data->suppressAlg!=1){ + if(data->suppressAlg!=0 && data->suppressAlg!=1){ printf("suppressAlg must either be 0 or 1!\n"); return(NULL); } @@ -421,7 +428,9 @@ UserData allocateUserData(FILE *input){ data->output=fopen("output.dat","w"); data->globalOutput=fopen("globalOutput.dat","w"); data->gridOutput=fopen("grid.dat","w"); - data->timescaleOutput=fopen("timeScale.dat","w") ; + data->timescaleOutput=fopen("timeScale.dat","w") ; + data->rxnROPOutput=fopen("rxnROP.dat","w"); + data->spROPOutput=fopen("spROP.dat","w"); //data->ratesOutput=fopen("rates.dat","w"); data->innerMassFractions = new double [data->nsp]; diff --git a/UserData.h b/UserData.h index 00379d6..ff713c4 100644 --- a/UserData.h +++ b/UserData.h @@ -209,6 +209,9 @@ typedef struct UserDataTag{ int JJRG; double deltaT; + FILE* rxnROPOutput; + FILE* spROPOutput; + } *UserData; UserData allocateUserData(FILE *input); diff --git a/main.cpp b/main.cpp index 820597b..ab0e15f 100644 --- a/main.cpp +++ b/main.cpp @@ -214,6 +214,7 @@ int main(){ int move=0; int kcur=0; int RGCOUNT=0; + size_t ii=0; if(data->adaptiveGrid){ dxMin=data->grid->leastMove; xOld=maxCurvPosition(ydata, data->nt, data->nvar, @@ -222,9 +223,6 @@ int main(){ // data->nvar, data->grid->x, data->npts); } while (tNow<=finalTime && R(1)>100e-9) { - - - t1=tNow; /*Floor small value to zero*/ @@ -257,7 +255,6 @@ int main(){ if(data->adaptiveGrid==1 && data->moveGrid==1){ x=maxCurvPosition(ydata, data->nt, data->nvar, data->grid->x, data->npts); - //x=isothermPosition(ydata, data->isotherm, data->nt, // data->nvar, data->grid->x, data->npts); @@ -397,8 +394,14 @@ int main(){ // printSpaceTimeRates(tNow, ydot, data); //} } - - + + /*Get and Print Rxns Rate of Progress and Specie Rate of Production data*/ + /*Following code snippet will be executed only once*/ + if(ii==0 && maxT >=(data->initialTemperature+data->deltaT)){ + getReactions(data,&y,data->rxnROPOutput); + getSpecies(data,&y,data->spROPOutput); + ii++; + } // getTimescale(data,&y); // if(count%data->nSaves==0){ // printTimescaleOutput(tNow,&y, data->timescaleOutput,data); diff --git a/residue.cpp b/residue.cpp index bdc654b..724caec 100644 --- a/residue.cpp +++ b/residue.cpp @@ -80,7 +80,6 @@ void REGRID(double* ydata,double* ydotdata,UserData data){ } INTERPO(ydata,ydotdata,nvar,nPts,XNEW,XOLD); - } @@ -148,7 +147,7 @@ double maxTemperature(const double* y,const size_t nt,const size_t nvar ,size_t } return(maxT); } - +//Index here is 1-based int maxTemperatureIndex(const double* y,const size_t nt,const size_t nvar ,size_t nPts){ double maxT = 0.0e0; double TempT = 0.0e0; @@ -671,7 +670,6 @@ int initializePsiGrid(double* ydata, double* psidata, UserData data){ int setInitialCondition(N_Vector* y, N_Vector* ydot, UserData data){ - double* ydata; double* ydotdata; double* psidata; @@ -1108,9 +1106,9 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ /*Update specific heat:Cpl for liquid proprane&n-heptane&dodecane*/ /*Based on the existiong data,a linear expression is used*/ // Cpl= 12.10476*T(1) - 746.60143; // Unit:J/(kg*K), Need to be improved later - Cpl[0] = 2.423*T(1)+1661.074; //Unit:J/(kg*K),range:1atm,propane - Cpl[1] = 3.336*T(1)+1289.5; //Unit:J/(kg*K),range:1atm,n-Heptane - //Cpl[1] = 3.815*T(1)+1081.8; //Unit:J/(kg*K),range:1atm,n-Dodecane + //Cpl[0] = 2.423*T(1)+1661.074; //Unit:J/(kg*K),range:1atm,propane + Cpl[0] = 3.336*T(1)+1289.5; //Unit:J/(kg*K),range:1atm,n-Heptane + Cpl[1] = 3.815*T(1)+1081.8; //Unit:J/(kg*K),range:1atm,n-Dodecane double Cp_l = 0.0; for(size_t i=0;i<=1;i++){ Cp_l=Cp_l+Cpl[i]*data->dropMassFrac[i]; @@ -1137,9 +1135,9 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ /*Following section is related to the property of liquid n-heptane and n-dodecane (mainly density rho)*/ /*Density of these two species should be temperature dependent*/ - data->dropDens[0] = -1.046*T(1)+823.794; //Unit:kg/m^3,density of propane @1atm - data->dropDens[1] = -0.858*T(1)+933.854; //Unit:kg/m^3,density of n-Heptane @1atm - //data->dropDens[1] = -0.782*T(1)+979.643; //Unit:kg/m^3,density of n-Dodecane @1atm + //data->dropDens[0] = -1.046*T(1)+823.794; //Unit:kg/m^3,density of propane @1atm + data->dropDens[0] = -0.858*T(1)+933.854; //Unit:kg/m^3,density of n-Heptane @1atm + data->dropDens[1] = -0.782*T(1)+979.643; //Unit:kg/m^3,density of n-Dodecane @1atm //rhol = data->dropRho; //Unit:kg/m^3 rhol = 0.0; for(size_t i=0;i<=1;i++){ @@ -1149,7 +1147,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ aream= calc_area(R(1),&m); - + /*******************************************************************/ /*Calculate values at j=2's m and mhalf*****************************/ @@ -1209,7 +1207,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ ///Yres(1,k)=innerMassFractionsData[k-1]- /// Y(1,k)- /// (YVmhalf(k)*areamhalf)/Mdot(1); - //Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf; + //} //else{ // //Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1]; @@ -1793,11 +1791,8 @@ void getTimescale(UserData data, N_Vector* y){ for(k=1;k<= nsp;k++){ data->time_scale(i,k) = concentra(k)/(wdot_mass(k)+ 1.00e-16) ; } - - } - } @@ -1836,7 +1831,6 @@ void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data) fprintf(output, "\n"); } fprintf(output, "\n"); - } void floorSmallValue(UserData data, N_Vector* y) @@ -1915,3 +1909,75 @@ void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv) } +void getReactions(UserData data,N_Vector* y,FILE* output){ + double Tmax; + double* ydata; + //double deltaT = 400.0; + //int i = 0; + size_t nRxns; + int index; + nRxns = data->gas->nReactions(); + //DEBUG + printf("Total Number of Rxns:%zu \n",nRxns); + double fwdROP[nRxns],revROP[nRxns],netROP[nRxns]; + std::string* rxnArr = new std::string[nRxns]; + /*DEBUG*/ + printf("Memory Allocation for Rxns Array Succeed!\n"); + ydata = N_VGetArrayPointer_OpenMP(*y); + Tmax = maxTemperature(ydata,data->nt,data->nvar,data->npts); + //get the rxns' equation array + for(size_t ii=0;iigas->reactionString(ii); + } + if(Tmax>=(data->initialTemperature+data->deltaT)){ + index = maxTemperatureIndex(ydata,data->nt,data->nvar,data->npts); + setGas(data,ydata,index); + /*Get forward/reverse/net rate of progress of rxns*/ + data->gas->getRevRatesOfProgress(revROP); + data->gas->getFwdRatesOfProgress(fwdROP); + data->gas->getNetRatesOfProgress(netROP); + for(size_t j=0 ; jnsp],revROP[data->nsp],netROP[data->nsp]; + std::string* spArr = new std::string[data->nsp]; + /*DEBUG*/ + printf("Memory Allocation for Species Array Succeed!\n"); + ydata = N_VGetArrayPointer_OpenMP(*y); + Tmax = maxTemperature(ydata,data->nt,data->nvar,data->npts); + //get the species name array + for(size_t ii=0;iinsp;ii++){ + spArr[ii]= data->gas->speciesName(ii); + } + + if(Tmax>=(data->initialTemperature+data->deltaT)){ + //i = i + 1 ; + index = maxTemperatureIndex(ydata,data->nt,data->nvar,data->npts); + setGas(data,ydata,index); + /*Get forward/reverse/net rate of progress of rxns*/ + data->gas->getDestructionRates(revROP); + data->gas->getCreationRates(fwdROP); + data->gas->getNetProductionRates(netROP); + /*Print data to the pertinent output file*/ + for(size_t j=0 ; jnsp; j++){ + fprintf(output,"%15s\t",spArr[j]); + fprintf(output,"%15.9e\t%15.9e\t%15.9e\t\n",fwdROP[j],revROP[j],netROP[j]); + } + fprintf(output, "\n"); + fclose(output); + } + delete[] spArr; +} diff --git a/residue.h b/residue.h index 349766c..eb29392 100644 --- a/residue.h +++ b/residue.h @@ -111,3 +111,6 @@ void printTimescaleHeader(UserData data); void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data); void floorSmallValue(UserData data, N_Vector* y); void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv); + +void getReactions(UserData data,N_Vector* y,FILE* output); +void getSpecies(UserData data,N_Vector* y,FILE* output);