diff --git a/DropletCombustion11 b/DropletCombustion11 new file mode 100755 index 0000000..2bcebce Binary files /dev/null and b/DropletCombustion11 differ diff --git a/DropletCombustionTest1 b/DropletCombustionTest1 new file mode 100755 index 0000000..1054c58 Binary files /dev/null and b/DropletCombustionTest1 differ diff --git a/DropletCombustionTest3 b/DropletCombustionTest3 new file mode 100755 index 0000000..fb8a1d2 Binary files /dev/null and b/DropletCombustionTest3 differ diff --git a/DropletCombustionTest4 b/DropletCombustionTest4 new file mode 100755 index 0000000..a5ff4fb Binary files /dev/null and b/DropletCombustionTest4 differ diff --git a/DropletCombustionTest5 b/DropletCombustionTest5 new file mode 100755 index 0000000..188d8f4 Binary files /dev/null and b/DropletCombustionTest5 differ diff --git a/Makefile b/Makefile index 19e235c..4b9a584 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 =DropletCombustion9 +EXE =DropletCombustionTest6 #DESTDIR =~/bin DESTDIR = ../example diff --git a/UserData.cpp b/UserData.cpp index c0aa148..5102224 100644 --- a/UserData.cpp +++ b/UserData.cpp @@ -380,6 +380,12 @@ UserData allocateUserData(FILE *input){ ier=reGrid(data->grid, data->grid->position); if(ier==-1)return(NULL); + +// /**************** TEST THE data->grid->xOld *******************/ +// double* ptr = data->grid->xOld ; +// printf("allocateUserData function is called,Start print the first 5 elements of the xOld array : \n"); +// printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr[0],ptr[1],ptr[2],ptr[3],ptr[4]); + data->npts=data->grid->nPts; data->neq=data->nvar*data->npts; } diff --git a/gridRoutines.cpp b/gridRoutines.cpp index 489e03d..6f6107b 100644 --- a/gridRoutines.cpp +++ b/gridRoutines.cpp @@ -1,5 +1,6 @@ #include "gridRoutines.h" #include + inline double l(const double* x, const double* a, const double* w, @@ -8,10 +9,40 @@ inline double l(const double* x, if(*refineLeft==0){ return(tanh(-*a*(*x+*w*100.0e0))); }else{ - return(tanh(-*a*(*x-*w*(*fac)))); + double l ; + l = tanh(-*a*(*x-*w*(*fac))) ; + + if(l>=0){ + return l*10.0; + }else{ + //return(tanh(-*a*(*x-*w*(*fac)))); + return l; + } } } +//inline double l(const double* x, +// const double* a, +// const double* c, +// const double* w, +// const double* fac, +// const int* refineLeft){ +// if(*refineLeft==0){ +// return(tanh(-*a*(*x+*w*100.0e0))); +// }else{ +// double l ; +// //l = tanh(-*a*(*x-*w*(*fac))) ; +// l = tanh(-*a*(*x-*c)); +// +// if(l>=0){ +// return l*10.0; +// }else{ +// //return(tanh(-*a*(*x-*w*(*fac)))); +// return l; +// } +// } +//} + inline double r(const double* x, const double* a, const double* w, @@ -55,6 +86,7 @@ inline double rho(const double* x, return(((2.0e0+f(x,a,c,w) +g(x,a,c,w) +l(x,a,w,leftFac,refineLeft) +// +l(x,a,c,w,leftFac,refineLeft) +r(x,a,w,rightFac,refineRight))*0.5e0) *(*mag-1.0e0)+1.0e0); } @@ -122,16 +154,33 @@ void fillGrid(const size_t* basePts, double sum=0.0e0; double err=0.0e0; double fix=0.0e0; + double arr[*nPts-1] ; + size_t halfboundII = 0; for(size_t j=0;j<*nPts-1;j++){ sum+=dxp[j]; + arr[j]=sum; } + + for(size_t j=0;j<*nPts;j++){ + if(arr[j] > 0.5e0){ + halfboundII = j; + break; + } + } + err=1.0e0-sum; printf("sum before correction: %15.6e\n",sum); printf("err before correction: %15.6e\n",err); - fix=err/((double)(*nPts)); + + //fix=err/((double)(*nPts)); + fix = err/((double)(*nPts-1-(halfboundII+1))); + sum=0.0e0; for(size_t j=0;j<*nPts-1;j++){ - dxp[j]+=fix; + // dxp[j]+=fix; + if(j>halfboundII){ + dxp[j]=dxp[j]+fix ; + } sum+=dxp[j]; } err=1.0e0-sum; diff --git a/gridRoutines.h b/gridRoutines.h index 06231ac..9352fd7 100644 --- a/gridRoutines.h +++ b/gridRoutines.h @@ -33,6 +33,13 @@ inline double l(const double* x, const double* fac, const int* refineLeft); +//inline double l(const double* x, +// const double* a, +// const double* c, +// const double* w, +// const double* fac, +// const int* refineLeft); + inline double r(const double* x, const double* a, const double* w, diff --git a/main.cpp b/main.cpp index 868bbb8..d97460a 100644 --- a/main.cpp +++ b/main.cpp @@ -38,6 +38,12 @@ int main(){ fclose(input); data->clockStart=get_wall_time(); + +// /**************** TEST THE xOld *******************/ +// double* ptr1 = data->grid->xOld ; +// printf("After allocateUserData in main.cpp,Start print the first 5 elements of the xOld array : \n"); +// printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr1[0],ptr1[1],ptr1[2],ptr1[3],ptr1[4]); + if(data==NULL){ printf("check input file!\n"); freeUserData(data); @@ -53,6 +59,13 @@ int main(){ y=ydot=id=res=atolv=constraints=NULL; ier=allocateSolution(data->neq,data->nThreads,&y,&ydot,&id,&res,&atolv,&constraints); ier=setAlgebraicVariables(&id,data); + + +// /**************** TEST THE xOld *******************/ +// double* ptr2 = data->grid->xOld ; +// printf("Before setInitialCondition in main.cpp,Start print the first 5 elements of the xOld array : \n"); +// printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr2[0],ptr2[1],ptr2[2],ptr2[3],ptr2[4]); + ier=setInitialCondition(&y,&ydot,data); if(ier==-1)return(-1); tNow=data->tNow; @@ -154,10 +167,16 @@ int main(){ if(!data->dryRun){ printf("Calculating Initial Conditions:\n"); printf("Cross your fingers...\n"); - ier = IDACalcIC(mem, IDA_YA_YDP_INIT, 1e-05*finalTime); + + ier = IDACalcIC(mem, IDA_YA_YDP_INIT, 1e-5*finalTime); + //If at first IDACalcIC doesn't succeed, try, try, try again: for (int i = 0; i < 10; i++) { ier = IDACalcIC(mem, IDA_YA_YDP_INIT, (1e-01+pow(10,i)*finalTime)); + + /************* Print the #of iterations ************/ + //printf("This the %dth try of calculating initial conditions. \n",i); + if(ier==0){ break; } @@ -179,8 +198,8 @@ int main(){ printSpaceTimeOutput(tNow, &y, data->output, data); printSpaceTimeOutput(tNow, &y, data->gridOutput, data); - getTimescale(data,&y) ; - printTimescaleOutput(tNow, &y, data->timescaleOutput,data); +// getTimescale(data,&y) ; +// printTimescaleOutput(tNow, &y, data->timescaleOutput,data); if(!data->dryRun){ @@ -196,13 +215,15 @@ int main(){ int kcur=0; if(data->adaptiveGrid){ dxMin=data->grid->leastMove; - //xOld=maxCurvPosition(ydata, data->nt, data->nvar, - // data->grid->x, data->npts); - xOld=isothermPosition(ydata, data->isotherm, data->nt, - data->nvar, data->grid->x, data->npts); + xOld=maxCurvPosition(ydata, data->nt, data->nvar, + data->grid->x, data->npts); + //xOld=isothermPosition(ydata, data->isotherm, data->nt, + // data->nvar, data->grid->x, data->npts); } while (tNow<=finalTime && R(1)>100e-9) { + + t1=tNow; /*Floor small value to zero*/ @@ -229,15 +250,20 @@ int main(){ dt=tNow-t1; ier = IDAGetCurrentOrder(mem, &kcur); + /******** Print the max Temperature *********/ + double maxT = 0.00; + maxT = maxTemperature(ydata,data->nt,data->nvar,data->npts); + printf("Maximum temperature is : %.3f[K] \n",maxT); + if(data->adaptiveGrid==1 && data->moveGrid==1){ - //x=maxCurvPosition(ydata, data->nt, data->nvar, - // data->grid->x, data->npts); + 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); + // data->nvar, data->grid->x, data->npts); - x=maxGradPosition(ydata, data->nt, data->nvar, - data->grid->x, data->npts); + //x=maxGradPosition(ydata, data->nt, data->nvar, + // data->grid->x, data->npts); dx=x-xOld; if(dx*dxMin>0.0e0){ @@ -248,6 +274,11 @@ int main(){ //if(fabs(dx)>=dxMin && x+(double)(move)*0.5e0*dxMin<=1.0e0){ dxRatio=fabs(dx/dxMin); + + /************ Print xOld, x,dx,dxMin,dxRatio ******************/ + printf("xOld = %.6f, x = %.6f,",xOld,x); + printf("dx = %.6f, dxMin = %.6f, dxRatio = %.3f\n",dx,dxMin,dxRatio); + if(dxRatio>=1.0e0 && dxRatio<=2.0e0){ printf("Regridding!\n"); @@ -319,14 +350,14 @@ int main(){ } - getTimescale(data,&y); - if(count%data->nSaves==0){ - printTimescaleOutput(tNow,&y, data->timescaleOutput,data); - //printSpaceTimeOutput(tNow, &y, data->output, data); - //if(data->writeRates){ - // printSpaceTimeRates(tNow, ydot, data); - //} - } +// getTimescale(data,&y); +// 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){ diff --git a/readme.md b/readme.md index 9d00db8..d4f408f 100644 --- a/readme.md +++ b/readme.md @@ -3,3 +3,12 @@ Dropletcombustion6 is the version where the Antoine parameters are those of wate Version 7 change the hardcode some tolerance after ignition:Max T> 1800K version 8: the x in the main.cpp is changed to the maxGradPosition version 9: isothermPosition function in residue.cpp is revised for droplet combustion +version 10: both xOld and x is the position of isoTherm ; +version 11: both xOld and x is the position of maximum curvature; +DropletCombustionTest: test the variables in the program +DropletCombustionTest1: Multiple the l function by 10.0 +DropletCombustionTest2: Multiple the l function by 100.0 and remove print timeScale; +DropletCombustionTest3: Multiple the l function by 10.0 and print the maxTemp and add criteria for regrid +DropletCombustionTest4: dxRatio to be 1.0e-2 and delta_T to be 150 [K] +DropletCombustionTest5: dxRatio to be 1.0 and print x and xOld temperature, modified the fillGrid function; + diff --git a/residue.cpp b/residue.cpp index c87a563..f4c4944 100644 --- a/residue.cpp +++ b/residue.cpp @@ -9,6 +9,115 @@ #include #include "timing.hpp" +double maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double *x ,size_t nPts){ + double maxT = 0.0e0; + double TempT = 0.0e0; + double pos = 0.0e0; + for (size_t i=1;i<=nPts;i++){ + TempT = y[(i-1)*nvar+nt] ; + if(TempT > maxT){ + maxT = TempT ; + pos = x[i-1]; + } + } + return(pos); +} + +double maxTemperature(const double* y,const size_t nt,const size_t nvar ,size_t nPts){ + double maxT = 0.0e0; + double TempT = 0.0e0; + //int index = 0 ; + for (size_t i=1;i<=nPts;i++){ + TempT = y[(i-1)*nvar+nt] ; + if(TempT > maxT){ + maxT = TempT ; + } + } + return(maxT); +} + +int maxTemperatureIndex(const double* y,const size_t nt,const size_t nvar ,size_t nPts){ + double maxT = 0.0e0; + double TempT = 0.0e0; + int index = 0 ; + for (size_t i=1;i<=nPts;i++){ + TempT = y[(i-1)*nvar+nt] ; + if(TempT > maxT){ + maxT = TempT ; + index = i; + } + } + return(index); +} + +double maxCurvPositionR(const double* y, const size_t nt, + const size_t nvar, const size_t nr, size_t nPts){ + double maxCurvT=0.0e0; + double gradTp=0.0e0; + double gradTm=0.0e0; + double curvT=0.0e0; + double dx=0.0e0; + double dr=0.0e0; + double pos=0.0e0; + size_t j,jm,jp; + size_t r,rp,rm; + for (size_t i = 1; i =maxCurvT) { + maxCurvT=curvT; + pos=y[r]; + } + } + return(pos); +} + +int maxCurvIndexR(const double* y, const size_t nt, + const size_t nvar, const size_t nr, size_t nPts){ + double maxCurvT=0.0e0; + double gradTp=0.0e0; + double gradTm=0.0e0; + double curvT=0.0e0; + double dx=0.0e0; + double dr=0.0e0; + int pos=0; + size_t j,jm,jp; + size_t r,rm,rp; + for (size_t i = 1; i =maxCurvT) { + maxCurvT=curvT; + pos=i; + } + } + return(pos); +} + double maxGradPosition(const double* y, const size_t nt, const size_t nvar, const double* x, size_t nPts){ double maxGradT=0.0e0; @@ -31,7 +140,7 @@ int maxGradIndex(const double* y, const size_t nt, const size_t nvar, const double* x, size_t nPts){ double maxGradT=0.0e0; double gradT=0.0e0; - int pos=0.0e0; + int pos=0; size_t j,jm; for (size_t i = 1; i domainLength)/((double)(data->npts)-1.0e0); double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0); - for (size_t i = 1; i <=data->npts; i++) { + + //Initialize the R(i),T(i)(data->initialTemperature),Y(i,k),P(i) + for (size_t i = 1; i <=data->npts; i++) { if(data->metric==0){ R(i)=Rd+(double)((i-1)*dR); }else{ @@ -497,6 +608,18 @@ int setInitialCondition(N_Vector* y, P(i)=data->initialPressure*Cantera::OneAtm; } R(data->npts)=data->domainLength+data->Rd; + +// /********** test R(i) and the volumn between grids *****************/ +// printf("Print the first 4 R(i)s and volumn between them: \n") ; +// printf("Grids: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e \n",R(1),R(2),R(3),R(4)); +// double v1,v2,v3,v4,v5; +// v1 =pow(R(2),1+data->metric) - pow(R(1),1+data->metric); +// v2 =pow(R(3),1+data->metric) - pow(R(2),1+data->metric); +// v3 =pow(R(4),1+data->metric) - pow(R(3),1+data->metric); +// v4 =pow(R(5),1+data->metric) - pow(R(4),1+data->metric); +// v5 =pow(R(6),1+data->metric) - pow(R(5),1+data->metric); +// printf("Volumn: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e,5th:%2.6e\n",v1,v2,v3,v4,v5) ; + double Tmax; double Tmin=data->initialTemperature; @@ -554,8 +677,21 @@ int setInitialCondition(N_Vector* y, data->grid->position=0.0e0; double x=data->grid->position+data->gridOffset*data->grid->leastMove; printf("New grid center:%15.6e\n",x); - ier=reGrid(data->grid, x); + +// /**************** TEST THE data->grid->xOld *******************/ +// double* ptr3 = data->grid->xOld ; +// printf("SetInitialCondition function is called,before reGrid,Start print the first 5 elements of the xOld array : \n"); +// printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr3[0],ptr3[1],ptr3[2],ptr3[3],ptr3[4]); + + ier=reGrid(data->grid, x); if(ier==-1)return(-1); + +// /**************** TEST THE data->grid->xOld *******************/ +// double* ptr = data->grid->xOld ; +// printf("SetInitialCondition function is called,after reGrid,Start print the first 5 elements of the xOld array : \n"); +// printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr[0],ptr[1],ptr[2],ptr[3],ptr[4]); + + updateSolution(ydata, ydotdata, data->nvar, data->grid->xOld,data->grid->x,data->npts); storeGrid(data->grid->x,data->grid->xOld,data->npts); @@ -662,16 +798,24 @@ int setInitialCondition(N_Vector* y, } } - //Set grid to location of maximum curvature: + //Set grid to location of maximum curvature calculated by r instead of x: if(data->adaptiveGrid){ - data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar, - data->grid->x, data->npts); + //data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar, + // data->grid->x, data->npts); + int maxCurvII = 0; + maxCurvII = maxCurvIndexR(ydata,data->nt,data->nvar,data->nr,data->npts); + + data->grid->position = data->grid->x[maxCurvII]; ier=reGrid(data->grid, data->grid->position); updateSolution(ydata, ydotdata, data->nvar, data->grid->xOld,data->grid->x,data->npts); storeGrid(data->grid->x,data->grid->xOld,data->npts); + + /******** Test the maxCurvPosition and related variables *******/ + printf("The maxCurvPositition(based on r) is : %.6f,Temperature is : %.3f \n",data->grid->position,T(maxCurvII+1)); } + /*Ensure consistent boundary conditions*/ //T(1)=T(2); //for (size_t k = 1; k <=data->nsp; k++) { diff --git a/residue.h b/residue.h index bd7d448..128733c 100644 --- a/residue.h +++ b/residue.h @@ -19,6 +19,17 @@ #include "UserData.h" +double maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double* x ,size_t nPts); + +double maxTemperature(const double* y,const size_t nt, const size_t nvar, size_t nPts); + +int maxTemperatureIndex(const double* y,const size_t nt,const size_t nvar ,size_t nPts); + +double maxCurvPositionR(const double* y, const size_t nt, + const size_t nvar, const size_t nr, size_t nPts); + +int maxCurvIndexR(const double* y, const size_t nt, + const size_t nvar, const size_t nr, size_t nPts); double maxGradPosition(const double* y, const size_t nt, const size_t nvar, const double* x, size_t nPts);