diff --git a/DropletCombustionTest7 b/DropletCombustionTest7 new file mode 100755 index 0000000..8534194 Binary files /dev/null and b/DropletCombustionTest7 differ diff --git a/Makefile b/Makefile index 4b9a584..dd67f3d 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 =DropletCombustionTest6 +EXE =DropletCombustionTest7-binary #DESTDIR =~/bin DESTDIR = ../example diff --git a/UserData.cpp b/UserData.cpp index 5102224..85b5de6 100644 --- a/UserData.cpp +++ b/UserData.cpp @@ -339,20 +339,29 @@ UserData allocateUserData(FILE *input){ } } - ier=parseNumber(input, "dropletComposition", MAXBUFLEN, data->dropSpec); - if(ier==-1){ - printf("Enter composition of droplet!\n"); - return(NULL); - } + //ier=parseNumber(input, "dropletComposition", MAXBUFLEN, data->dropSpec); + //if(ier==-1){ + // printf("Enter composition of droplet!\n"); + // return(NULL); + //} + ier=parseDrop(input,"dropletComposition",data->dropSpec,data->dropMole,MAXBUFLEN); + ier=parseDrop(input,"dropletDensity",data->dropSpec,data->dropDens,MAXBUFLEN); - ier=parseNumber (input, "nThreads" , MAXBUFLEN, &data->nThreads); + + ier=parseNumber (input, "nThreads", MAXBUFLEN, &data->nThreads); if(data->nThreads<0 ){ printf("nThreads must be greater than 0!\n"); return(NULL); } + - data->nr=0; + ier=parseNumber(input, "PCAD", MAXBUFLEN, &data->PCAD); + ier=parseNumber(input,"RGTC", MAXBUFLEN, &data->RGTC); + ier=parseNumber(input,"JJRG", MAXBUFLEN, &data->JJRG); + ier=parseNumber(input,"deltaT", MAXBUFLEN, &data->deltaT); + + data->nr=0; //data->np=data->nr+1; data->nt=data->nr+1; data->ny=data->nt+1; @@ -361,9 +370,28 @@ UserData allocateUserData(FILE *input){ data->nvar=data->nsp+4; //assign no: of variables (R,T,P,Mdot,nsp species) + + double MW[2]; + for(int i=0;i<=1;i++){ + int speciesIndex = data->gas->speciesIndex(data->dropSpec[i]); + double weight = data->gas->molecularWeight(speciesIndex); + MW[i]= weight; + } + + double massSum = 0.0; + for(int i=0;i<=1;i++){ + massSum = massSum + data->dropMole[i] * MW[i]; + } + data->dropRho = 0.0; + for(int i=0;i<=1;i++){ + data->dropMassFrac[i] = data->dropMole[i]*MW[i]/massSum; + data->dropRho = data->dropRho + data->dropMassFrac[i]*data->dropDens[i]; + } + printf("Density of droplet is: %.3f [kg/m^3]\n",data->dropRho); //Mass of droplet //data->massDrop=1.0/3.0*pow(data->Rd,3)*997.0; //TODO: The density of the droplet should be a user input - data->massDrop=1.0/3.0*pow(data->Rd,3)*684.00; //TODO: The density of the droplet(n-heptane) should be a user input + //data->massDrop=1.0/3.0*pow(data->Rd,3)*684.00; //TODO: The density of the droplet(n-heptane) should be a user input + data->massDrop=1.0/3.0*pow(data->Rd,3)*data->dropRho; if(!data->adaptiveGrid){ data->uniformGrid = new double [data->npts]; @@ -454,5 +482,9 @@ void setSaneDefaults(UserData data){ data->flameTime[0]=0.0e0; data->flameTime[1]=0.0e0; data->nTimeSteps=0; + data->PCAD=0.75; + data->RGTC=1.0; + data->JJRG=0; + data->deltaT=200.0; } diff --git a/UserData.h b/UserData.h index da595f3..00379d6 100644 --- a/UserData.h +++ b/UserData.h @@ -5,6 +5,7 @@ #endif #include "gridRoutines.h" +#include #ifndef USER_DEF #define USER_DEF @@ -16,7 +17,15 @@ typedef struct UserDataTag{ * species.*/ Cantera::Transport* trmix; /* Droplet species composition */ - char dropSpec[MAXBUFLEN]; + //char dropSpec[MAXBUFLEN]; + char dropSpec[2][10]; + /* Droplet species mole fractions */ + double dropMole[2]; + /* Droplet species density at given initialTemperature*/ + double dropDens[2]; + /* Droplet species mass fractions*/ + double dropMassFrac[2]; + double dropRho; /*Length of the domain (in meters):*/ double domainLength; /*Initial Droplet Radius (in meters)*/ @@ -54,7 +63,8 @@ typedef struct UserDataTag{ size_t k_OH; size_t k_HO2; /*Species index of droplet composition*/ - size_t k_drop; + /*Index starts with 1 instead of 0*/ + size_t k_drop[2]; /*User-defined mass flux (kg/m^2/s):*/ double mdot; /*Flag to solve isobaric/isochoric problem;*/ @@ -193,6 +203,12 @@ typedef struct UserDataTag{ //double* MW; FILE* timescaleOutput; + /*Following parameters are used for REGRID function*/ + double PCAD; + double RGTC; + int JJRG; + double deltaT; + } *UserData; UserData allocateUserData(FILE *input); diff --git a/macros.h b/macros.h index 3c827b8..bcfbc8c 100644 --- a/macros.h +++ b/macros.h @@ -15,6 +15,9 @@ * macros here. Also, we define macros to ease referencing various variables in * the sundials nvector. */ +#define WORK(i) WORK[i-1] +#define XNEW(i) XNEW[i-1] + #define psi(i) psidata[i-1] #define T(i) ydata[((i-1)*data->nvar)+data->nt] diff --git a/main.cpp b/main.cpp index d97460a..820597b 100644 --- a/main.cpp +++ b/main.cpp @@ -213,6 +213,7 @@ int main(){ double dxRatio=dx/dxMin; int move=0; int kcur=0; + int RGCOUNT=0; if(data->adaptiveGrid){ dxMin=data->grid->leastMove; xOld=maxCurvPosition(ydata, data->nt, data->nvar, @@ -229,8 +230,6 @@ int main(){ /*Floor small value to zero*/ floorSmallValue(data, &y); - /*reset the tolerance after ignition*/ - resetTolerance(data,&y,&atolv); if(data->quasiSteady){ ier = IDASolve(mem, finalTime, &tNow, y, ydot, IDA_ONE_STEP); @@ -338,6 +337,56 @@ int main(){ } } + + /*reset the tolerance after ignition*/ + //resetTolerance(data,&y,&atolv); + + /*regrid and update the solution based on R,re-initialize the problem*/ + /*For the time being,we only allow TORC to REGRID once for each run*/ + + if(data->JJRG ==1 && (maxT >= data->initialTemperature+data->deltaT)){ + if(RGCOUNT<1){ + RGCOUNT = RGCOUNT +1; + REGRID(ydata,ydotdata,data); + initializePsiGrid(ydata,data->uniformGrid,data); + printf("REGRID Complete!Restarting Problem at %15.6e s\n",tNow); + ier = IDAReInit(mem,tNow,y,ydot); + + if(check_flag(&ier,"IDAReInit",1)){ + freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data); + return(-1); + } + + ier= IDASetInitStep(mem,1e-01*dt); + //ier= IDASetInitStep(mem,0.0); + printf("Reinitialized!Calculating Initial Conditions:\n"); + printf("Cross your fingers...\n"); + ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e-01*dt); + if(check_flag(&ier, "IDACalcIC", 1)){ + ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e+01*dt); + } + //Every once in a while, for reasons + //that befuddle this mathematically + //lesser author, IDACalcIC fails. Here, + //I desperately try to make it converge + //again by sampling increasingly larger + //time-steps: + for (int i = 0; i < 10; i++) { + ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+(1e-01+pow(10,i)*dt)); + if(ier==0){ + break; + } + } + //Failure :( Back to the drawing board: + if(check_flag(&ier, "IDACalcIC", 1)){ + freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data); + return(-1); + } + printf("Initial (Consistent) Conditions Calculated!\n"); + printf("------------------------------------------\n\n"); + + } + } /*Floor small value to zero*/ floorSmallValue(data, &y); diff --git a/parse.cpp b/parse.cpp index f872ca5..d2e6e9f 100644 --- a/parse.cpp +++ b/parse.cpp @@ -18,3 +18,62 @@ void getFromString (const char* buf, char* n){ sscanf(buf,"%s",n); printf("%s\n",n); } + +/*Extract droplet species and mole fractions*/ +int parseDrop(FILE* input, const char* keyword,char dropSpec[][10],double dropMole[],const size_t bufLen){ + char buf[bufLen]; + char buf1[bufLen]; + char comment[1]; + char *ret; + while (fgets(buf,bufLen, input)!=NULL){ + comment[0]=buf[0]; + if(strncmp(comment,"#",1)==0){ + } + else{ + strcpy(buf1,buf); + ret=strtok(buf,"="); + //DEBUG + //printf("Current KEYWORD in input: %20s \n",ret); + if(strcmp(ret,keyword)==0){ + char* modifiedFuel = NULL; + char* equalSign = strstr(buf1,"="); + if(equalSign!= NULL){ + modifiedFuel = new char [strlen(equalSign)+1]; + strcpy(modifiedFuel,equalSign+1); + //DEBUG + //printf("modifiedFuel:%20s \n",modifiedFuel); + char* token = strtok(modifiedFuel,","); + int index = 0 ; + char* list[2]; + while(token!= NULL){ + //DEBUG + //printf("TOKEN :%20s \n",token); + list[index] = token; + token = strtok(NULL,","); + index++; + } + //for (int i=0;i<2;i++){ + // printf("%20s",list[i]); + //} + for(int i=0;i<2;i++){ + char* name= strtok(list[i],":"); + char* value = strtok(NULL,":"); + //DEBUG + // printf("Name:%10s,Value:%10s \n",name,value); + strcpy(dropSpec[i],name); + dropMole[i]=std::stod(value); + // printf("In the dropArray,Name:%10s,Value:%.3f\n",dropSpec[i],dropMole[i]); + } + delete[] modifiedFuel; + } + printf("%10s:%10s:%.3f,%10s:%.3f\n",keyword,dropSpec[0],dropMole[0],dropSpec[1],dropMole[1]); + //printf("IF statement is execuated. \n"); + rewind(input); + //delete[] modifiedFuel; + return(0); + } + } + } + rewind(input); + return(-1); +} diff --git a/parse.h b/parse.h index 45f3a2c..894bd12 100644 --- a/parse.h +++ b/parse.h @@ -3,6 +3,10 @@ #include //for strings #include //for printf,scanf #include //for atoi, atof +#include +#include +#include +#include #endif #ifndef PARSE_DEF @@ -23,6 +27,7 @@ template int parseArray(FILE* input, const char* keyword, const size_t bufLen, const size_t arrLen, T y[]); +int parseDrop(FILE* input, const char* keyword,char dropSpec[][10],double dropMole[],const size_t bufLen); #include "parse.hpp" diff --git a/parse.hpp b/parse.hpp index 634db96..6c10a9a 100644 --- a/parse.hpp +++ b/parse.hpp @@ -73,4 +73,3 @@ int parseArray(FILE* input, const char* keyword, const size_t bufLen, rewind(input); return(-1); } - diff --git a/readme.md b/readme.md index d4f408f..8f798f4 100644 --- a/readme.md +++ b/readme.md @@ -11,4 +11,6 @@ DropletCombustionTest2: Multiple the l function by 100.0 and remove print timeSc 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; +DropletCombustionTest7: JJRG has been implemented; + diff --git a/residue.cpp b/residue.cpp index f4c4944..bdc654b 100644 --- a/residue.cpp +++ b/residue.cpp @@ -9,6 +9,119 @@ #include #include "timing.hpp" +void REGRID(double* ydata,double* ydotdata,UserData data){ + size_t nPts = data->npts; + size_t nvar = data->nvar; + double WORK[nPts],XNEW[nPts],XOLD[nPts]; + double P,R,R0,R1,R2,TV1,TV2,XLEN,B1,B2; + double DX,DDX; + double DETA,ETAJ,DEL; + size_t ISTART; + P = data->PCAD; + R = data->RGTC; + R0 =1.0 - P; + R1 =P*R/(R+1.0); + R2 =P - R1 ; + TV1 = 0.0; + TV2 = 0.0; + + for(size_t i=2;i<=nPts;i++){ + TV1 = TV1 + fabs(T(i)-T(i-1)); + } + + for(size_t i=2;i<=nPts-1;i++){ + TV2 = TV2 + fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) ); + } + + XLEN = R(nPts) - R(1); + B1 = R1* XLEN/ ( (1.0-P)*TV1); + B2 = R2* XLEN/ ( (1.0-P)*TV2); + + //Compute partial sum of weight function + WORK(1) =0.0; + DX = 0.0; + for(size_t i=2;i<=nPts-1;i++){ + DX = R(i) - R(i-1) ; + WORK(i) = DX + B1*fabs(T(i)-T(i-1)) + WORK(i-1)+B2*fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) ); + } + + DDX = R(nPts) - R(nPts-1); + WORK(nPts) = WORK(nPts-1) +DDX ; + + for(size_t i=2;i<=nPts;i++){ + WORK(i) = WORK(i)/WORK(nPts); + } + + XNEW(1)=R(1); + XNEW(nPts)=R(nPts); + ISTART =2; + //double DETA; + DETA = 1.0/(nPts-1) ; + + //Fill new grid XNEW[nPts],duplicate from PREMIX REGRID SUBROUTINE + for(size_t j=2;j<=nPts-1;j++){ + ETAJ= (j-1)*DETA; + + for(size_t i=ISTART;i<=nPts;i++){ + if(ETAJ <= WORK(i)){ + DEL = (ETAJ-WORK(i-1))/(WORK(i)-WORK(i-1)) ; + XNEW(j)=R(i-1)+(R(i)-R(i-1))*DEL; + break; + }else{ + ISTART =i; + } + + } + } + + //Interpolate solution based on XNEW[]&XOLD[] + for(size_t i=1;i<=nPts;i++){ + XOLD[i-1]=R(i); + } + + INTERPO(ydata,ydotdata,nvar,nPts,XNEW,XOLD); + +} + + +void INTERPO(double* y,double* ydot,const size_t nvar,size_t nPts,const double XNEW[], const double XOLD[]){ + double ytemp[nPts],ydottemp[nPts] ; + + gsl_interp_accel* acc; + gsl_spline* spline; + acc = gsl_interp_accel_alloc(); + spline = gsl_spline_alloc(gsl_interp_steffen,nPts); + + gsl_interp_accel* accdot; + gsl_spline* splinedot; + accdot = gsl_interp_accel_alloc(); + splinedot = gsl_spline_alloc(gsl_interp_steffen,nPts); + + for(size_t j=0;jk_bath=BathGasIndex(data); data->k_oxidizer=oxidizerIndex(data); data->k_OH=OHIndex(data); data->k_HO2=HO2Index(data); - data->k_drop=specIndex(data,data->dropSpec); + //use char* pointer to get the address + for(int i=0;i<2;i++){ + specPtr[i] = data->dropSpec[i]; + } + data->k_drop[0]=specIndex(data,specPtr[0]); + data->k_drop[1]=specIndex(data,specPtr[1]); + printf("Oxidizer index: %lu\n",data->k_oxidizer); printf("Bath gas index:%lu\n",data->k_bath); - printf("Droplet species index:%lu\n",data->k_drop); + printf("Droplet species index:%lu,%lu\n",data->k_drop[0],data->k_drop[1]); for (size_t i = 1; i <=data->npts; i++) { /*Algebraic variables: indicated by ZERO.*/ Rid(i)=ZERO; @@ -361,7 +482,8 @@ int setAlgebraicVariables(N_Vector* id, UserData data){ } Mdotid(1)=ZERO; Rid(1)=ONE; - Yid(1,data->k_drop)=ZERO; + Yid(1,data->k_drop[0])=ZERO; + Yid(1,data->k_drop[1])=ZERO; Tid(data->npts)=ZERO; if(data->constantPressure){ Pid(data->npts)=ONE; @@ -569,8 +691,9 @@ int setInitialCondition(N_Vector* y, //Mass of droplet //massDrop=1.0/3.0*Rd*Rd*Rd*997.0; //TODO: The density of the droplet should be a user input - massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input - if(data->adaptiveGrid){ + //massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input + massDrop = 1.0/3.0*Rd*Rd*Rd*data->dropRho; + if(data->adaptiveGrid){ psidata = data->grid->xOld; } else{ @@ -903,7 +1026,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ size_t npts=data->npts; size_t nsp=data->nsp; size_t k_bath = data->k_bath; - size_t k_drop = data->k_drop; + size_t k_drop[2] = {data->k_drop[0],data->k_drop[1]}; ydata = N_VGetArrayPointer_OpenMP(y); ydotdata= N_VGetArrayPointer_OpenMP(ydot); @@ -952,7 +1075,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ double rhomhalf, rhom, lambdamhalf, YVmhalf[nsp], rho, rhophalf, lambdaphalf, YVphalf[nsp], - Cpb, Cvb, Cp[nsp], Cpl, dHvl, rhol, wdot[nsp], enthalpy[nsp], energy[nsp], + Cpb, Cvb, Cp[nsp], Cpl[2], dHvl[2], rhol, wdot[nsp], enthalpy[nsp], energy[nsp], tranTerm, diffTerm, srcTerm, advTerm, area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq; @@ -982,24 +1105,46 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ //Cpl= 5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2) //J/kg/K for liquid water // - 9.62429538e+01*T(1) + 1.27131046e+04; - /*Update specific heat:Cpl for liquid n-heptane*/ + /*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 = 2242.871642; //Unit:J/(kg*K),range:25~520K,Ginnings&Furukawa,NIST + 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 + double Cp_l = 0.0; + for(size_t i=0;i<=1;i++){ + Cp_l=Cp_l+Cpl[i]*data->dropMassFrac[i]; + } //dHvl=2.260e6; //J/kg heat of vaporization of water //double Tr = T(1)/647.1; //Reduced Temperature: Droplet Temperature divided by critical temperature of water //dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01; //J/kg latent heat of vaporization of water - - double Tr = T(1)/540.2; //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST - dHvl = 5.366e5*exp(-0.2831*Tr) * pow((1-Tr),0.2831); //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST + double Tr1 = T(1)/369.8; //Reduced Temperature; + dHvl[0] = 6.32716e5*exp(-0.0208*Tr1)*pow((1-Tr1),0.3766); //Unit:J/kg,Latent Heat of Vaporizaiton of Liquid n-propane,NIST + + double Tr2 = T(1)/540.2; //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST + dHvl[1] = 5.366e5*exp(-0.2831*Tr2) * pow((1-Tr2),0.2831); //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST + + double dHv_l = 0.0; + for(size_t i=0;i<=1;i++){ + dHv_l=dHv_l+dHvl[i]*data->dropMassFrac[i]; + } + /*Following section is related to the property of water*/ //rhol = 997.0; //massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input - /*Following section is related to the property of liquid n-heptane (mainly density rho)*/ - rhol = 684.0; //Unit:kg/m^3 + /*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 + //rhol = data->dropRho; //Unit:kg/m^3 + rhol = 0.0; + for(size_t i=0;i<=1;i++){ + rhol = rhol + data->dropMassFrac[i]*data->dropDens[i]; + } massDrop = 1.0/3.0*R(1)*R(1)*R(1)*rhol; //Unit:kg, this is the mass of liquid droplet aream= calc_area(R(1),&m); @@ -1038,18 +1183,24 @@ 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))) ; //unit:Pa,FOR N-HEPTANE + double p_i[2]={0.0,0.0} ; + + p_i[0] = 1.0e5*pow(10,4.53678-(1149.36/(T(1)+24.906)))*data->dropMassFrac[0] ; //unit:Pa,For N-PROPANE + p_i[1] = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199)))*data->dropMassFrac[1] ; //unit:Pa,FOR N-HEPTANE //p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0))); //Unit:Pa,FOR WATER - Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) - / P(1) / data->gas->meanMolecularWeight(); - - sum=sum+Y(1,k_drop); + //Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) + // / P(1) / data->gas->meanMolecularWeight(); + for(size_t i=0;i<=1;i++){ + Yres(1,k_drop[i])=Y(1,k_drop[i])-p_i[i]*data->gas->molecularWeight(k_drop[i]-1)/(P(1)*data->gas->meanMolecularWeight()); + } - Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass + + sum=sum+Y(1,k_drop[0])+Y(1,k_drop[1]); + //Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass + Mdotres(1)= Mdot(1) - ((YVmhalf(k_drop[0])+YVmhalf(k_drop[1])) *areamhalf) / (1-Y(1,k_drop[0])-Y(1,k_drop[1]) ); //Evaporating mass for (k = 1; k <=nsp; k++) { - if(k!=k_bath && k!=k_drop){ + if(k!=k_bath && k!=k_drop[0] && k!=k_drop[1]){ //TODO: May need to include chemical source term Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf; //Yres(1,k)=YVmhalf(k)*areamhalf; @@ -1076,14 +1227,19 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){ if(data->dirichletInner){ //Tres(1)=Tdot(1); //Tres(1)=T(1)-data->innerTemperature; - Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHvl - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop)); + /*Following code should be revised in the future*/ + //Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHv_l - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop)); }else{ //Tres(1)=Tdot(1) - // (rhomhalf*areamhalfsq*lambdamhalf*(T(2)-T(1))/((psi(2)-psi(1))*mass) - Mdot(1)*dHvl) // / (massDrop * Cpl); + + //Tres(1)=Tdot(1) - + // (areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl) + // / (massDrop * Cpl); Tres(1)=Tdot(1) - - (areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl) - / (massDrop * Cpl); + (areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHv_l) + / (massDrop * Cp_l); @@ -1720,25 +1876,25 @@ void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv) ydata = N_VGetArrayPointer_OpenMP(*y); atolvdata = N_VGetArrayPointer_OpenMP(*atolv); double relTol,radTol,tempTol,presTol,massFracTol,bathGasTol,mdotTol; - double maxT=0.00, ignT=1800; + double maxT=0.00, deltaT=400.0; relTol = 1.0e-03 ; radTol = 1.0e-05 ; - tempTol = 1.0e-01 ; - presTol = 1.0e00 ; - massFracTol = 1.0e-08; - bathGasTol = 1.0e-04 ; - mdotTol = 1.0e-09 ; + tempTol = 1.0e-03 ; + presTol = 1.0e-3 ; + massFracTol = 1.0e-06; + bathGasTol = 1.0e-05 ; + mdotTol = 1.0e-10 ; - /*Get the maximum Tempture*/ + /*Get the maximum Temperature*/ for (size_t i =1;i <= data->npts;i++){ if(T(i) > maxT){ maxT = T(i); } } - /*reset the tolerance when maxT > ignT*/ - if(maxT >= ignT){ + /*reset the tolerance when maxT > initialTemperature +deltaT*/ + if(maxT >= (data->initialTemperature + deltaT)){ data->relativeTolerance = relTol; for(size_t i =1; i<=data->npts; i++){ diff --git a/residue.h b/residue.h index 128733c..349766c 100644 --- a/residue.h +++ b/residue.h @@ -7,6 +7,8 @@ #ifndef PRINT_DEF #define PRINT_DEF #include //for strings +#include +#include #include //for printf,scanf #include //for atoi, atof #endif @@ -19,6 +21,10 @@ #include "UserData.h" +void REGRID(double* ydata,double* ydotdata,UserData data); + +void INTERPO(double* y,double* ydot,const size_t nvar,size_t nPts,const double XNEW[], const double XOLD[]); + 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);