Liquid Phase Property of n-Heptane/n-Dodecane updated

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
 

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<int>   (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];

UserData.h

@@ -209,6 +209,9 @@ typedef struct UserDataTag{
     int JJRG;
     double deltaT;
 
+    FILE* rxnROPOutput;
+    FILE* spROPOutput;
+
 
 } *UserData;
 UserData allocateUserData(FILE *input);

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);

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;ii<nRxns;ii++){
+        rxnArr[ii]= data->gas->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 ; j<nRxns; j++){
+            fprintf(output,"%30s\t",rxnArr[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[] rxnArr;
+}
+
+void getSpecies(UserData data,N_Vector* y,FILE* output){
+    double Tmax;
+    double* ydata; 
+    //double deltaT = 400.0;
+    //int i = 0;
+    int index;
+    double fwdROP[data->nsp],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;ii<data->nsp;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 ; j<data->nsp; 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;
+}

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);