diff --git a/Makefile b/Makefile index cf0bc60..8f2c047 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 =DropletCombustion4 +EXE =DropletCombustion5 #DESTDIR =~/bin DESTDIR = ../example diff --git a/main.cpp b/main.cpp index 3c69992..a43151c 100644 --- a/main.cpp +++ b/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){ diff --git a/readme.md b/readme.md new file mode 100644 index 0000000..e95074d --- /dev/null +++ b/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 diff --git a/residue.cpp b/residue.cpp index 7c1f65f..8a691da 100644 --- a/residue.cpp +++ b/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; + } + } diff --git a/residue.h b/residue.h index 4b168b7..c532ec4 100644 --- a/residue.h +++ b/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); \ No newline at end of file +void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data); +void floorSmallValue(UserData data, N_Vector* y);