function floorSmallValue added to main.cpp

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		=DropletCombustion4
+EXE		=DropletCombustion5
 #DESTDIR		=~/bin
 DESTDIR		= ../example
 

main.cpp

@@ -204,6 +204,9 @@ int main(){
 		while (tNow<=finalTime && R(1)>100e-9) {
 
 			t1=tNow;
+            
+            /*Floor small value to zero*/
+            floorSmallValue(data, &y); 
 
 			if(data->quasiSteady){
   				ier = IDASolve(mem, finalTime, &tNow, y, ydot, IDA_ONE_STEP);
@@ -302,6 +305,9 @@ int main(){
 				}
 			}
 
+            /*Floor small value to zero*/
+            floorSmallValue(data, &y);
+
 			if(count%data->nSaves==0 && !data->writeEveryRegrid){
 				printSpaceTimeOutput(tNow, &y, data->output, data);
 				//if(data->writeRates){

readme.md

@@ -0,0 +1,3 @@
+Dropletcombustion5 is the version which floor the small mass fraction(smaller than massfractionTolerance) to 0.00 ; 
+Dropletcombustion6 is the version where the Antoine parameters are those of water ;
+Version 7 change the hardcode some tolerance after ignition:Max T> 1800K

residue.cpp

@@ -880,7 +880,8 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
     /*Following using the Antoine Parameters to calculate the pressure of volatile component(n-heptane)*/
     //Antoine parameter from NIST website
     double p_i=0.0 ; 
-    p_i = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199))) ;
+    p_i = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199))) ;    //FOR N-HEPTANE 
+    //p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)));     //FOR WATER
     Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) 
 				    / P(1) / data->gas->meanMolecularWeight();
 	
@@ -1512,7 +1513,8 @@ void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data)
 
 	for(size_t i=1 ; i<= data->npts; i++){
 		fprintf(output, "%15.9e\t%15.9e\t",t,R(i));
-
+		fprintf(output, "%15.9e\t",T(i));
+        
 		for(size_t k=1; k<=data->nsp;k++){
 			fprintf(output, "%15.9e\t",data->time_scale(i,k));
 		}
@@ -1520,5 +1522,34 @@ void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data)
 	}
 	fprintf(output, "\n");
 
+}
+
+void floorSmallValue(UserData data, N_Vector* y)
+{
+    double* ydata;
+    ydata = N_VGetArrayPointer_OpenMP(*y);
+    //double sum = 0.00; 
+    size_t k_bath = data->k_bath;
+   
+	/*Floor small values to zero*/
+	for (size_t i = 1; i <=data->npts; i++) {
+		for (size_t k = 1; k <=data->nsp; k++) {
+			if(fabs(Y(i,k))<=data->massFractionTolerance){
+				Y(i,k)=0.0e0;
+			}
+		}
+	}
+    
+    /*Dump the error to the bath gas*/
+	for (size_t i = 1; i <=data->npts; i++) {
+        double sum = 0.00;
+		for (size_t k = 1; k <=data->nsp; k++) {
+			if(k!=k_bath){
+				sum = sum + Y(i,k);
+			}
+		}
+        Y(i,k_bath) = ONE-sum;
+	}
+    
 
 }

residue.h

@@ -91,4 +91,5 @@ 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);
+void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data);
+void floorSmallValue(UserData data, N_Vector* y);