Characteristic time scales

UserData.cpp

@@ -53,6 +53,10 @@ void freeUserData(UserData data){
   			fclose(data->globalOutput);
 			printf("Global Output File Cleared from Memory!\n");
 		}
+		if(data->timescaleOutput!=NULL){
+  			fclose(data->timescaleOutput);
+			printf("Characteristic Timescale Output File Cleared from Memory!\n");
+		}
 
 	}
   	free(data);             /* Free the user data */
@@ -383,9 +387,14 @@ 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->ratesOutput=fopen("rates.dat","w");
 	
 	data->innerMassFractions = new double [data->nsp];
+	//data->wdot_mole = new double [data->nsp] ; 
+	//data->wdot_mass = new double [data->nsp] ; 
+	data->time_scale = new double [(data->npts) * (data->nsp)] ; 
+	//data->MW = new double [data->nsp] ; 
 
 	return(data);
 }

UserData.h

@@ -184,6 +184,16 @@ typedef struct UserDataTag{
 	double flamePosition[2];
 	double flameTime[2];
 	size_t nTimeSteps;
+
+	/*Following arrays are used to compute the characteristic time scale of
+	species*/
+	//double* wdot_mole ; 
+	//double* wdot_mass ; 
+	double* time_scale ; 
+	//double* MW;
+	FILE* timescaleOutput;
+
+
 } *UserData;
 UserData allocateUserData(FILE *input);
 void setSaneDefaults(UserData data);

macros.h

@@ -66,4 +66,12 @@
 #define constraintsY(i,k)   constraintsdata[((i-1)*data->nvar)+data->ny+k-1]
 #define constraintsR(i)   constraintsdata[((i-1)*data->nvar)+data->nr]
 
+/*Following marcos are defined to calculate the characteristic time-scale */
+#define wdot_mole(i)    wdot_mole[i-1]
+#define wdot_mass(i)    wdot_mass[i-1]
+#define MW(i)         MW[i-1]
+#define time_scale(i,k)    time_scale[(i-1)*data->nsp+k-1]
+#define concentra(i)     concentra[i-1]
+
+
 #endif

main.cpp

@@ -175,8 +175,11 @@ int main(){
 
 	printSpaceTimeHeader(data);
 	printGlobalHeader(data);
+	printTimescaleHeader(data);
 	printSpaceTimeOutput(tNow, &y, data->output, data);
 	printSpaceTimeOutput(tNow, &y, data->gridOutput, data);
+	printTimescaleOutput(tNow, &y, data->timescaleOutput,data);
+	
 
 	if(!data->dryRun){
 		count=0;
@@ -304,6 +307,14 @@ int main(){
 				//}
 			}
 
+			if(count%data->nSaves==0 ){
+				printTimescaleOutput(tNow,&y, data->timescaleOutput,data);
+				//printSpaceTimeOutput(tNow, &y, data->output, data);
+				//if(data->writeRates){
+				//	printSpaceTimeRates(tNow, ydot, data);
+				//}
+			}
+
 			/*Print global variables only if time-step is of high order!*/
 			if(data->nTimeSteps==0){
 				data->flamePosition[0]=0.0e0;

residue.cpp

@@ -1444,3 +1444,80 @@ void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)
 ////		    Ydot(npts-1,k)=Ydot(npts,k);
 ////	}
 ////}
+
+
+/*Following functions are added to derive the characteristic time scale of species*/
+void getTimescale(UserData data, N_Vector* y){
+	size_t i, k, nsp,npts ; 
+	nsp = data->nsp ; 
+	npts = data->npts ; 
+	double rho, wdot_mole[nsp], wdot_mass[nsp],MW[nsp],concentra[nsp]; 
+	//double time_scale[npts*nsp] ;   
+
+	double* ydata ; 
+	ydata = N_VGetArrayPointer_OpenMP(*y); //access the data stored in N_Vector* y 
+	
+	for(i=1; i<= npts;i++){
+		setGas(data,ydata,i); //set the gas state at each grid point
+		rho = data->gas->density(); //get the averaged density at each grid point, Unit:kg/m^3
+		data->gas->getNetProductionRates(wdot_mole) ; //Unit:kmol/(m^3 * s)
+		
+		for(k=1; k<= nsp; k++){
+			MW(k) = data->gas->molecularWeight(k-1) ; //Unit:kg/kmol
+		}
+
+		for(k=1; k<= nsp; k++){
+			wdot_mass(k) = wdot_mole(k) * MW(k) ; //Unit:kg/(m^3 * s)
+		}
+
+		for(k=1; k<= nsp; k++){
+			concentra(k) = Y(i,k) * rho ;  //Unit:kg/m^3
+		}
+
+		for(k=1;k<= nsp;k++){
+			data->time_scale(i,k) = concentra(k)/(wdot_mass(k)+ 1.00e-16) ; 
+		}
+
+
+	}
+
+}
+
+
+void printTimescaleHeader(UserData data)
+{
+	fprintf((data->timescaleOutput), "%15s\t","#1");
+	for (size_t k = 1; k <=data->nsp+1; k++) {
+		fprintf((data->timescaleOutput), "%15lu\t",k+1);
+	}
+	fprintf((data->timescaleOutput), "%15lu\n",data->nsp+3);
+
+
+	fprintf((data->timescaleOutput), "%15s\t%15s\t%15s\t","#time","radius","Temp(K)");
+	//fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
+	for (size_t k = 1; k <=(data->nsp); k++) {
+		fprintf((data->timescaleOutput), "%15s\t",data->gas->speciesName(k-1).c_str());
+	}
+	//fprintf((data->output), "%15s\t","Pressure(Pa)");
+	//fprintf((data->timescaleOutput), "%15s\n",data->gas->speciesName(data->nsp-1).c_str());
+	//fprintf((data->output), "%15s\n","Mdot (kg/s)");
+	fprintf((data->timescaleOutput), "\n");
+}
+
+void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data)
+{
+	double* ydata; 
+	ydata = N_VGetArrayPointer_OpenMP(*y); 
+
+	for(size_t i=1 ; i<= data->npts; i++){
+		fprintf(output, "%15.9e\t%15.9e\t",t,R(i));
+
+		for(size_t k=1; k<=data->nsp;k++){
+			fprintf(output, "%15.9e\t",data->time_scale(i,k));
+		}
+		fprintf(output, "\n");
+	}
+	fprintf(output, "\n");
+
+
+}

residue.h

@@ -88,3 +88,7 @@ void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data);
 
 void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data);
 void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data);
+
+void getTimescale(UserData data, N_Vector* y);
+void printTimescaleHeader(UserData data);
+void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data);