weiye
/
Droplet_LTORC
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
Droplet Lagrangian Transient One-dimensional Reacting Code Implementation of both liquid and gas phase governing equations.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
6 Commits
1 Branch
1.6MB
Tree: d28fc466b5
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'd28fc466b5'
${ noResults }
Droplet_LTORC/residue.cpp

1618 lines
48KB
Raw Blame History 

  1. #ifndef GSL_DEF

  2. #define GSL_DEF

  3. #include <gsl/gsl_math.h>

  4. #include <gsl/gsl_spline.h>

  5. #endif

  6. #include "residue.h"

  7. #include "macros.h"

  8. #include <cmath>

  9. #include <stdio.h>

  10. #include "timing.hpp"

  11. 

  12. double maxGradPosition(const double* y, const size_t nt, 

  13. 	               const size_t nvar, const double* x, size_t nPts){

  14. 	double maxGradT=0.0e0;

  15. 	double gradT=0.0e0;

  16. 	double pos=0.0e0;

  17. 	size_t j,jm;

  18. 	for (size_t i = 1; i <nPts; i++) {

  19. 		j=i*nvar+nt;

  20. 		jm=(i-1)*nvar+nt;

  21. 		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));

  22. 		if (gradT>=maxGradT) {

  23. 			maxGradT=gradT;

  24. 			pos=x[i];

  25. 		}

  26. 	}

  27. 	return(pos);

  28. }

  29. 

  30. int maxGradIndex(const double* y, const size_t nt, 

  31. 	         const size_t nvar, const double* x, size_t nPts){

  32. 	double maxGradT=0.0e0;

  33. 	double gradT=0.0e0;

  34. 	int pos=0.0e0;

  35. 	size_t j,jm;

  36. 	for (size_t i = 1; i <nPts; i++) {

  37. 		j=i*nvar+nt;

  38. 		jm=(i-1)*nvar+nt;

  39. 		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));

  40. 		if (gradT>=maxGradT) {

  41. 			maxGradT=gradT;

  42. 			pos=i;

  43. 		}

  44. 	}

  45. 	return(pos);

  46. }

  47. 

  48. double maxCurvPosition(const double* y, const size_t nt, 

  49. 	               const size_t nvar, const double* x, size_t nPts){

  50. 	double maxCurvT=0.0e0;

  51. 	double gradTp=0.0e0;

  52. 	double gradTm=0.0e0;

  53. 	double curvT=0.0e0;

  54. 	double dx=0.0e0;

  55. 	double pos=0.0e0;

  56. 	size_t j,jm,jp;

  57. 	for (size_t i = 1; i <nPts-1; i++) {

  58. 		j=i*nvar+nt;

  59. 		jm=(i-1)*nvar+nt;

  60. 		jp=(i+1)*nvar+nt;

  61. 		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));

  62. 		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));

  63. 		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));

  64. 		curvT=(gradTp-gradTm)/dx;

  65. 		if (curvT>=maxCurvT) {

  66. 			maxCurvT=curvT;

  67. 			pos=x[i];

  68. 		}

  69. 	}

  70. 	return(pos);

  71. }

  72. 

  73. int maxCurvIndex(const double* y, const size_t nt, 

  74. 	               const size_t nvar, const double* x, size_t nPts){

  75. 	double maxCurvT=0.0e0;

  76. 	double gradTp=0.0e0;

  77. 	double gradTm=0.0e0;

  78. 	double curvT=0.0e0;

  79. 	double dx=0.0e0;

  80. 	int pos=0;

  81. 	size_t j,jm,jp;

  82. 	for (size_t i = 1; i <nPts-1; i++) {

  83. 		j=i*nvar+nt;

  84. 		jm=(i-1)*nvar+nt;

  85. 		jp=(i+1)*nvar+nt;

  86. 		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));

  87. 		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));

  88. 		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));

  89. 		curvT=(gradTp-gradTm)/dx;

  90. 		if (curvT>=maxCurvT) {

  91. 			maxCurvT=curvT;

  92. 			pos=i;

  93. 		}

  94. 	}

  95. 	return(pos);

  96. }

  97. 

  98. //double isothermPosition(const double* y, const double T, const size_t nt, 

  99. //	                const size_t nvar, const double* x, const size_t nPts){

  100. //	double pos=x[nPts-1];

  101. //	size_t j;

  102. //	for (size_t i = 1; i <nPts; i++) {

  103. //		j=i*nvar+nt;

  104. //		if (y[j]<=T) {

  105. //			pos=x[i];

  106. //			break;

  107. //		}

  108. //	}

  109. //	return(pos);

  110. //}

  111. 

  112. //******************* scan the temperature from left to right ******************//

  113. double isothermPosition(const double* y, const double T, const size_t nt, 

  114. 	                const size_t nvar, const double* x, const size_t nPts){

  115. 	double pos=x[0];

  116. 	size_t j;

  117. 	for (size_t i = 0; i <nPts; i++) {

  118. 		j=i*nvar+nt;

  119. 		if (y[j]>=T) {

  120. 			pos=x[i];

  121. 			break;

  122. 		}

  123. 	}

  124. 	return(pos);

  125. }

  126. 

  127. void updateSolution(double* y, double* ydot, const size_t nvar,

  128. 		    const double xOld[],const double xNew[],const size_t nPts){

  129. 

  130. 	double ytemp[nPts],ydottemp[nPts];

  131. 

  132. 	gsl_interp_accel* acc;

  133. 	gsl_spline* spline;

  134. 	acc = gsl_interp_accel_alloc();

  135. 	spline = gsl_spline_alloc(gsl_interp_steffen, nPts);

  136. 

  137. 	gsl_interp_accel* accdot;

  138. 	gsl_spline* splinedot;

  139. 	accdot    = gsl_interp_accel_alloc();

  140. 	splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);

  141. 

  142. 	for (size_t j = 0; j < nvar; j++) {

  143. 

  144. 		for (size_t i = 0; i < nPts; i++) {

  145. 			ytemp[i]=y[j+i*nvar];

  146. 			ydottemp[i]=ydot[j+i*nvar];

  147. 		}

  148. 

  149. 		gsl_spline_init(spline,xOld,ytemp,nPts);

  150. 		gsl_spline_init(splinedot,xOld,ydottemp,nPts);

  151. 

  152. 		for (size_t i = 0; i < nPts; i++) {

  153. 

  154. 		      	y[j+i*nvar]=gsl_spline_eval(spline,xNew[i],acc);

  155. 		      	ydot[j+i*nvar]=gsl_spline_eval(splinedot,xNew[i],accdot);

  156. 		}

  157. 	}

  158. 

  159. 	//Exploring "fixing" boundary conditions:

  160. 	//for (size_t j = 1; j < nvar; j++) {

  161. 	//	//printf("%15.6e\t%15.6e\n", y[j],y[j+nvar]);

  162. 	//	y[j]=y[j+nvar];

  163. 	//	//y[j+(nPts-1)*nvar]=y[j+(nPts-2)*nvar];

  164. 	//	//ydot[j+nvar]=ydot[j];

  165. 	//}

  166. 	//y[0]=0.0e0;

  167. 	

  168. 	gsl_interp_accel_free(acc);

  169. 	gsl_spline_free(spline);

  170. 

  171. 	gsl_interp_accel_free(accdot);

  172. 	gsl_spline_free(splinedot);

  173. }

  174. 

  175. //Locate bath gas:

  176. size_t BathGasIndex(UserData data){

  177. 	size_t index=0;

  178. 	double max1;

  179. 	double max=data->gas->massFraction(0);

  180. 	for(size_t k=1;k<data->nsp;k++){

  181. 		max1=data->gas->massFraction(k);

  182. 		if(max1>=max){

  183. 			max=max1;

  184. 			index=k;

  185. 		}

  186. 	}

  187. 	return(index+1);

  188. }

  189. 

  190. //Locate Oxidizer:

  191. size_t oxidizerIndex(UserData data){

  192. 	size_t index=0;

  193. 	for(size_t k=1;k<data->nsp;k++){

  194. 		if(data->gas->speciesName(k-1)=="O2"){

  195. 			index=k;

  196. 		}

  197. 	}

  198. 	return(index);

  199. }

  200. 

  201. //Locate OH:

  202. size_t OHIndex(UserData data){

  203. 	size_t index=0;

  204. 	for(size_t k=1;k<data->nsp;k++){

  205. 		if(data->gas->speciesName(k-1)=="OH"){

  206. 			index=k;

  207. 		}

  208. 	}

  209. 	return(index);

  210. }

  211. 

  212. //Locate HO2:

  213. size_t HO2Index(UserData data){

  214. 	size_t index=0;

  215. 	for(size_t k=1;k<data->nsp;k++){

  216. 		if(data->gas->speciesName(k-1)=="HO2"){

  217. 			index=k;

  218. 		}

  219. 	}

  220. 	return(index);

  221. }

  222. 

  223. //Locate species index:

  224. size_t specIndex(UserData data,char *specName){

  225. 	size_t index=0;

  226. 	for(size_t k=1;k<data->nsp;k++){

  227. 		if(data->gas->speciesName(k-1)==specName){

  228. 			index=k;

  229. 		}

  230. 	}

  231. 	return(index);

  232. }

  233. 

  234. int setAlgebraicVariables(N_Vector* id, UserData data){

  235. 	double *iddata;

  236.   	N_VConst(ONE, *id);

  237.   	iddata   = N_VGetArrayPointer_OpenMP(*id);

  238. 	data->k_bath=BathGasIndex(data);

  239. 	data->k_oxidizer=oxidizerIndex(data);

  240. 	data->k_OH=OHIndex(data);

  241. 	data->k_HO2=HO2Index(data);

  242. 	data->k_drop=specIndex(data,data->dropSpec);

  243. 	printf("Oxidizer index: %lu\n",data->k_oxidizer);

  244. 	printf("Bath gas index:%lu\n",data->k_bath);

  245. 	printf("Droplet species index:%lu\n",data->k_drop);

  246. 	for (size_t i = 1; i <=data->npts; i++) {

  247. 		/*Algebraic variables: indicated by ZERO.*/

  248. 		Rid(i)=ZERO;

  249. 		Pid(i)=ZERO;

  250. 		Yid(i,data->k_bath)=ZERO;

  251. 		Mdotid(i)=ZERO;

  252. 	}	

  253. 	Mdotid(1)=ZERO;

  254. 	Rid(1)=ONE;

  255. 	Yid(1,data->k_drop)=ZERO;

  256. 	Tid(data->npts)=ZERO;

  257. 	if(data->constantPressure){

  258. 		Pid(data->npts)=ONE;

  259. 	}else{

  260. 		Pid(data->npts)=ZERO;

  261. 		Rid(data->npts)=ONE;

  262. 	}

  263. 	if(data->dirichletInner){

  264. 		Tid(1)=ZERO; 

  265. 	}else{

  266. 		Tid(1)=ONE;

  267. 	}

  268. 	if(data->dirichletOuter){

  269. 		for (size_t k = 1; k <=data->nsp; k++) {

  270. 			if(k!=data->k_bath){

  271. 				Yid(data->npts,k)=ONE;

  272. 			}

  273. 			Yid(1,k)=ZERO;

  274. 			Tid(data->npts)=ONE;

  275. 		}

  276. 	}else{

  277. 		for (size_t k = 1; k <=data->nsp; k++){

  278. 			Yid(1,k)=ZERO;

  279. 			Yid(data->npts,k)=ZERO;

  280. 			Tid(data->npts)=ZERO;

  281. 		}

  282. 	}

  283. 	return(0);

  284. }

  285. 

  286. inline double calc_area(double x,int* i){

  287. 	switch (*i) {

  288. 		case 0:

  289. 			return(ONE);

  290. 		case 1:

  291. 			return(x);

  292. 		case 2:

  293. 			return(x*x);

  294. 		default:

  295. 			return(ONE);

  296. 	}

  297. }

  298. 

  299. void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts, double Rg){

  300. 

  301. 	FILE* output;output=fopen("test.dat","w");

  302. 

  303. 	size_t bufLen=10000;

  304. 	size_t nRows=0;

  305. 	size_t nColumns=nvar;

  306. 

  307. 	char buf[bufLen];

  308. 	char buf1[bufLen];

  309. 	char comment[1];

  310. 	char *ret;

  311. 

  312. 	while (fgets(buf,bufLen, input)!=NULL){

  313. 		comment[0]=buf[0];

  314. 		if(strncmp(comment,"#",1)!=0){

  315. 			nRows++;

  316. 		}

  317. 	}

  318. 	rewind(input);

  319. 

  320. 	printf("nRows: %ld\n", nRows);

  321. 	double y[nRows*nColumns];

  322. 

  323. 	size_t i=0;

  324. 	while (fgets(buf,bufLen, input)!=NULL){

  325. 		comment[0]=buf[0];

  326. 		if(strncmp(comment,"#",1)==0){

  327. 		}

  328. 		else{

  329. 			ret=strtok(buf,"\t");

  330. 			size_t j=0;

  331. 			y[i*nColumns+j]=(double)(atof(ret));

  332. 			j++;

  333. 			while(ret!=NULL){

  334. 				ret=strtok(NULL,"\t");

  335. 				if(j<nColumns){

  336. 					y[i*nColumns+j]=(double)(atof(ret));

  337. 				}

  338. 				j++;

  339. 			}

  340. 			i++;

  341. 		}

  342. 	}

  343. 

  344. 	for (i = 0; i < nRows; i++) {

  345. 		for (size_t j = 0; j < nColumns; j++) {

  346. 			fprintf(output, "%15.6e\t",y[i*nColumns+j]);

  347. 		}

  348. 		fprintf(output, "\n");

  349. 	}

  350. 	fclose(output);

  351. 

  352. 	//double xOld[nRows],xNew[nPts],ytemp[nPts];

  353. 	double xOld[nRows],xNew[nPts],ytemp[nRows];

  354. 

  355. 	for (size_t j = 0; j < nRows; j++) {

  356. 		xOld[j]=y[j*nColumns+nr];

  357. 	}

  358. 

  359. 	//double dx=(xOld[nRows-1] - xOld[0])/((double)(nPts)-1.0e0);

  360. 	//for (size_t j = 0; j < nPts-1; j++) {

  361. 	//	xNew[j]=(double)j*dx + xOld[0];

  362. 	//}

  363. 	//xNew[nPts-1]=xOld[nRows-1];

  364. 	double dX0;

  365. 	int NN = nPts-2;

  366. 	

  367. 	dX0 = (xOld[nRows-1]-xOld[0])*(pow(Rg,1.0/NN) - 1.0)/(pow(Rg,(NN+1.0)/NN) - 1.0);

  368. 	xNew[0] = xOld[0];

  369. 	for (size_t i = 1; i < nPts-1; i++) {

  370. 		xNew[i]=xNew[i-1]+dX0*pow(Rg,(i-1.0)/NN);

  371. 	}

  372. 	xNew[nPts-1] = xOld[nRows-1];

  373. 

  374. 	gsl_interp_accel* acc;

  375. 	gsl_spline* spline;

  376. 	acc = gsl_interp_accel_alloc();

  377. 	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);

  378. 

  379. 	for (size_t j = 0; j < nColumns; j++) {

  380. 

  381. 		for (size_t k = 0; k < nRows; k++) {

  382. 			ytemp[k]=y[j+k*nColumns];

  383. 		}

  384. 

  385. 		gsl_spline_init(spline,xOld,ytemp,nRows);

  386. 

  387. 		for (size_t k = 0; k < nPts; k++) {

  388. 		      	ydata[j+k*nColumns]=gsl_spline_eval(spline,xNew[k],acc);

  389. 		}

  390. 	}

  391. 	

  392. 	gsl_interp_accel_free(acc);

  393. 	gsl_spline_free(spline);

  394. 

  395. }

  396. 

  397. double systemMass(double* ydata,UserData data){

  398. 	double mass=0.0e0;

  399. 	double rho;

  400. 	for (size_t i = 2; i <=data->npts; i++) {

  401.   		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   

  402. 		rho=data->gas->density();                                           

  403. 		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&m);

  404. 		mass+=rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);

  405. 	}

  406. 	return(mass);

  407. }

  408. 

  409. 

  410. int initializePsiGrid(double* ydata, double* psidata, UserData data){

  411. 

  412. 	double rho,rhom;

  413. 	/*Create a psi grid that corresponds CONSISTENTLY to the spatial grid

  414. 	 * "R" created above. Note that the Lagrangian variable psi has  units

  415. 	 * of kg. */

  416. 	psi(1)=ZERO;

  417. 	for (size_t i = 2; i <=data->npts; i++) {

  418.   		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   

  419. 		rho=data->gas->density();                                           

  420.   		data->gas->setState_TPY(T(i-1), P(i-1), &Y(i-1,1));   

  421. 		rhom=data->gas->density();                                           

  422. 		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&data->metric);

  423. 		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);

  424. 		psi(i)=psi(i-1)+(R(i)-R(i-1))*(rho*calc_area(R(i),&data->metric)+rhom*calc_area(R(i-1),&data->metric))/TWO;

  425. 	}

  426. 

  427. 	/*The mass of the entire system is the value of psi at the last grid

  428. 	 * point. Normalize psi by this mass so that it varies from zero to

  429. 	 * one. This makes psi dimensionless. So the mass needs to be

  430. 	 * multiplied back in the approporiate places in the governing

  431. 	 * equations so that units match.*/

  432. 	data->mass=psi(data->npts);

  433. 	for (size_t i = 1; i <=data->npts; i++) {

  434. 		psi(i)=psi(i)/data->mass;

  435. 	}

  436. 	return(0);

  437. }

  438. 

  439. 

  440. int setInitialCondition(N_Vector* y, 

  441. 	    		N_Vector* ydot, 

  442.             		UserData data){

  443. 

  444. 	double* ydata;

  445. 	double* ydotdata;

  446. 	double* psidata;

  447. 	double* innerMassFractionsData, Rd, massDrop;

  448. 	double rhomhalf, lambdamhalf, YVmhalf[data->nsp];

  449. 	double f=ZERO;

  450. 	double g=ZERO;

  451. 

  452. 	double perturb,rho;

  453. 	double epsilon=ZERO;

  454. 	int m,ier;

  455.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  456.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  457. 	innerMassFractionsData =  data->innerMassFractions;

  458. 	Rd = data->Rd;

  459. 	massDrop = data->massDrop;

  460. 

  461. 	//Mass of droplet

  462. 	//massDrop=1.0/3.0*Rd*Rd*Rd*997.0; //TODO: The density of the droplet should be a user input

  463.     massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input 

  464. 	if(data->adaptiveGrid){

  465.   		psidata  = data->grid->xOld;

  466. 	}

  467. 	else{

  468.   		psidata  = data->uniformGrid;

  469. 	}

  470. 

  471. 	m=data->metric;

  472. 

  473. 	data->innerTemperature=data->initialTemperature;

  474. 	for (size_t k = 1; k <=data->nsp; k++) {

  475.   		innerMassFractionsData[k-1]=data->gas->massFraction(k-1);

  476. 	}

  477. 

  478. 	//Define Grid:

  479. 	double dR=(data->domainLength)/((double)(data->npts)-1.0e0);

  480. 	double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0);

  481. 	for (size_t i = 1; i <=data->npts; i++) {

  482. 		if(data->metric==0){

  483. 			R(i)=Rd+(double)((i-1)*dR);

  484. 		}else{

  485. 			if(i==1){

  486. 				R(i)=ZERO;

  487. 			}else if(i==2){

  488. 				R(i)=data->firstRadius*data->domainLength;

  489. 			}else{

  490. 				R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));

  491. 			}

  492. 		}

  493. 		T(i)=data->initialTemperature;

  494. 		for (size_t k = 1; k <=data->nsp; k++) {

  495.   	        	Y(i,k)=data->gas->massFraction(k-1);	//Indexing different in Cantera

  496. 		}

  497.   	        P(i)=data->initialPressure*Cantera::OneAtm;

  498. 	}	

  499. 	R(data->npts)=data->domainLength+data->Rd;

  500. 

  501. 	double Tmax;

  502. 	double Tmin=data->initialTemperature;

  503. 	double w=data->mixingWidth;

  504. 	double YN2=ZERO;

  505. 	double YO2=ZERO;

  506. 	double YFuel,YOxidizer,sum;

  507. 

  508. 	//if(data->problemType==0){

  509. 	//	data->gas->equilibrate("HP");

  510. 	//	data->maxTemperature=data->gas->temperature();

  511. 	//}

  512. 	//else if(data->problemType==1){

  513. 	//	/*Premixed Combustion: Equilibrium products comprise ignition

  514. 	//	 * kernel at t=0. The width of the kernel is "mixingWidth"

  515. 	//	 * shifted by "shift" from the center.*/

  516. 	//	data->gas->equilibrate("HP");

  517. 	//	Tmax=data->gas->temperature();

  518. 	//	for (size_t i = 1; i <=data->npts; i++) {

  519. 	//		g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x

  520. 	//		f=ONE-g;					//decreasing function of x

  521. 	//		T(i)=(Tmax-Tmin)*f+Tmin;

  522. 	//		for (size_t k = 1; k <=data->nsp; k++) {

  523. 	//			Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);

  524. 	//		}

  525. 	//	}

  526. 	//	if(data->dirichletOuter){

  527. 	//		T(data->npts)=data->wallTemperature;

  528. 	//	}

  529. 	//}

  530. 	//else if(data->problemType==2){

  531. 	//	FILE* input;

  532. 	//	if(input=fopen("initialCondition.dat","r")){

  533. 	//		readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);

  534. 	//		fclose(input);

  535. 	//	}

  536. 	//	else{

  537. 	//		printf("file initialCondition.dat not found!\n");

  538. 	//		return(-1);

  539. 	//	}

  540. 	//}

  541. 

  542. 	initializePsiGrid(ydata,psidata,data);

  543. 

  544. 	if(data->adaptiveGrid){

  545. 	//	if(data->problemType!=0){

  546. 	//		data->grid->position=maxGradPosition(ydata, data->nt, data->nvar,

  547. 	//						 data->grid->xOld, data->npts);

  548. 	//		//data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,

  549. 	//		//				 data->grid->xOld, data->npts);

  550. 	//	}

  551. 	//	else{

  552. 	//	}

  553. 		if(data->problemType!=3){

  554. 			data->grid->position=0.0e0;

  555. 			double x=data->grid->position+data->gridOffset*data->grid->leastMove;

  556. 			printf("New grid center:%15.6e\n",x);

  557. 			ier=reGrid(data->grid, x);

  558. 			if(ier==-1)return(-1);

  559. 			updateSolution(ydata, ydotdata, data->nvar,

  560. 			               data->grid->xOld,data->grid->x,data->npts);

  561. 			storeGrid(data->grid->x,data->grid->xOld,data->npts);

  562. 		}

  563. 	}

  564. 	else{

  565. 		double Rg = data->Rg, dpsi0;

  566. 		int NN = data->npts-2;

  567. 		double psiNew[data->npts];

  568. 		

  569. 		dpsi0 = (pow(Rg,1.0/NN) - 1.0)/(pow(Rg,(NN+1.0)/NN) - 1.0);

  570. 		psiNew[0] = 0;

  571. 		for (size_t i = 1; i < data->npts-1; i++) {

  572. 			psiNew[i]=psiNew[i-1]+dpsi0*pow(Rg,(i-1.0)/NN);

  573. 		}

  574. 		psiNew[data->npts-1] = 1.0;

  575. 		printf("Last point:%15.6e\n",psiNew[data->npts-1]);

  576. 		updateSolution(ydata, ydotdata, data->nvar,

  577. 		               data->uniformGrid,psiNew,data->npts);

  578. 		storeGrid(psiNew,data->uniformGrid,data->npts);

  579. 

  580. 		//double psiNew[data->npts];

  581. 		//double dpsi=1.0e0/((double)(data->npts)-1.0e0);

  582. 		//for (size_t i = 0; i < data->npts; i++) {

  583. 		//	psiNew[i]=(double)(i)*dpsi;

  584. 		//}

  585. 		//printf("Last point:%15.6e\n",psiNew[data->npts-1]);

  586. 		//updateSolution(ydata, ydotdata, data->nvar,

  587. 		//               data->uniformGrid,psiNew,data->npts);

  588. 		//storeGrid(psiNew,data->uniformGrid,data->npts);

  589. 	}

  590. 

  591. 	if(data->problemType==0){

  592. 		data->gas->equilibrate("HP");

  593. 		data->maxTemperature=data->gas->temperature();

  594. 	}

  595. 	else if(data->problemType==1){

  596. 		/*Premixed Combustion: Equilibrium products comprise ignition

  597. 		 * kernel at t=0. The width of the kernel is "mixingWidth"

  598. 		 * shifted by "shift" from the center.*/

  599. 		data->gas->equilibrate("HP");

  600. 		Tmax=data->gas->temperature();

  601. 		for (size_t i = 1; i <=data->npts; i++) {

  602. 			g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x

  603. 			f=ONE-g;					//decreasing function of x

  604. 			T(i)=(Tmax-Tmin)*f+Tmin;

  605. 			for (size_t k = 1; k <=data->nsp; k++) {

  606. 				Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);

  607. 			}

  608. 		}

  609. 		if(data->dirichletOuter){

  610. 			T(data->npts)=data->wallTemperature;

  611. 		}

  612. 	}

  613. 	else if(data->problemType==2){

  614. 		FILE* input;

  615. 		if(input=fopen("initialCondition.dat","r")){

  616. 			readInitialCondition(input, ydata, data->nvar, data->nr, data->npts, data->Rg);

  617. 			fclose(input);

  618. 		}

  619. 		else{

  620. 			printf("file initialCondition.dat not found!\n");

  621. 			return(-1);

  622. 		}

  623. 		initializePsiGrid(ydata,psidata,data);

  624. 	}

  625. 	else if(data->problemType==3){

  626. 		FILE* input;

  627. 		if(input=fopen("restart.bin","r")){

  628. 			readRestart(y, ydot, input, data);

  629. 			fclose(input);

  630. 			printf("Restart solution loaded!\n");

  631. 			printf("Problem starting at t=%15.6e\n",data->tNow);

  632. 			return(0);

  633. 		}

  634. 		else{

  635. 			printf("file restart.bin not found!\n");

  636. 			return(-1);

  637. 		}

  638. 	}

  639. 

  640. 	if(data->reflectProblem){

  641. 		double temp;

  642. 		int j=1;

  643. 		while (data->npts+1-2*j>=0) {

  644. 			temp=T(j);

  645. 			T(j)=T(data->npts+1-j);

  646. 			T(data->npts+1-j)=temp;

  647. 			for (size_t k = 1; k <=data->nsp; k++) {

  648. 				temp=Y(j,k);

  649. 				Y(j,k)=Y(data->npts+1-j,k);

  650. 				Y(data->npts+1-j,k)=temp;

  651. 			}

  652. 			j=j+1;

  653. 		}

  654. 	}

  655. 

  656. 	/*Floor small values to zero*/

  657. 	for (size_t i = 1; i <=data->npts; i++) {

  658. 		for (size_t k = 1; k <=data->nsp; k++) {

  659. 			if(fabs(Y(i,k))<=data->massFractionTolerance){

  660. 				Y(i,k)=0.0e0;

  661. 			}

  662. 		}

  663. 	}

  664. 

  665. 	//Set grid to location of maximum curvature:

  666. 	if(data->adaptiveGrid){

  667. 		data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,

  668. 				  		data->grid->x, data->npts);

  669. 		ier=reGrid(data->grid, data->grid->position);

  670. 		updateSolution(ydata, ydotdata, data->nvar,

  671. 		               data->grid->xOld,data->grid->x,data->npts);

  672. 		storeGrid(data->grid->x,data->grid->xOld,data->npts);

  673. 	}

  674. 

  675. 	/*Ensure consistent boundary conditions*/	

  676. 	//T(1)=T(2);

  677. 	//for (size_t k = 1; k <=data->nsp; k++) {

  678. 	//	Y(1,k)=Y(2,k);

  679. 	//	Y(data->npts,k)=Y(data->npts-1,k);

  680. 	//}

  681. 

  682. 	return(0);

  683. }

  684. 

  685. inline double Qdot(double* t, 

  686. 	    	   double* x, 

  687. 	    	   double* ignTime, 

  688. 	    	   double* kernelSize, 

  689. 	    	   double* maxQdot){

  690. 	double qdot;

  691. 	if(*x<=*kernelSize){

  692. 		if((*t)<=(*ignTime)){

  693. 			qdot=(*maxQdot);

  694. 		}

  695. 		else{

  696. 			qdot=0.0e0;

  697. 		}

  698. 	}else{

  699. 		qdot=0.0e0;

  700. 	}

  701. 	return(qdot);

  702. }

  703. 

  704. inline void setGas(UserData data, 

  705. 	    double *ydata, 

  706. 	    size_t gridPoint){

  707.   	data->gas->setTemperature(T(gridPoint)); 	        

  708.   	data->gas->setMassFractions_NoNorm(&Y(gridPoint,1)); 	

  709.   	data->gas->setPressure(P(gridPoint));

  710. }

  711. 

  712. void getTransport(UserData data, 

  713. 		  double *ydata, 

  714. 		  size_t gridPoint,

  715. 		  double *rho,

  716. 		  double *lambda,

  717. 		  double YV[]){

  718. 

  719. 	double YAvg[data->nsp],

  720. 		 XLeft[data->nsp],

  721. 		 XRight[data->nsp],

  722. 		 gradX[data->nsp];

  723. 

  724. 	setGas(data,ydata,gridPoint);

  725. 	data->gas->getMoleFractions(XLeft);                                     

  726. 	setGas(data,ydata,gridPoint+1);

  727. 	data->gas->getMoleFractions(XRight);                                     

  728.                                                                              

  729. 	for (size_t k = 1; k <=data->nsp; k++) {                                          

  730. 		YAvg(k)=HALF*(Y(gridPoint,k)+

  731. 			      Y(gridPoint+1,k));                                 

  732. 		gradX(k)=(XRight(k)-XLeft(k))/

  733. 			 (R(gridPoint+1)-R(gridPoint));                       

  734. 	}                                                                    

  735. 	double TAvg = HALF*(T(gridPoint)+T(gridPoint+1));

  736. 	double gradT=(T(gridPoint+1)-T(gridPoint))/

  737. 		       (R(gridPoint+1)-R(gridPoint));                                  

  738. 

  739. 

  740.   	data->gas->setTemperature(TAvg); 	        

  741.   	data->gas->setMassFractions_NoNorm(YAvg); 	

  742.   	data->gas->setPressure(P(gridPoint));

  743. 

  744. 	*rho=data->gas->density();                                       

  745. 	*lambda=data->trmix->thermalConductivity();                      

  746. 	data->trmix->getSpeciesFluxes(1,&gradT,data->nsp,                     

  747. 				      gradX,data->nsp,YV);                

  748. 	//setGas(data,ydata,gridPoint);

  749. }

  750. 

  751. int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){

  752. 

  753. 	/*Declare and fetch nvectors and user data:*/

  754. 

  755.   	double *ydata, *ydotdata, *resdata, *psidata, *innerMassFractionsData;

  756.   	

  757. 	UserData data;

  758. 	data = (UserData)user_data;

  759. 	size_t npts=data->npts;

  760. 	size_t nsp=data->nsp;

  761. 	size_t k_bath = data->k_bath;

  762. 	size_t k_drop = data->k_drop;

  763. 

  764. 	ydata 	= N_VGetArrayPointer_OpenMP(y); 

  765. 	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 

  766. 	resdata = N_VGetArrayPointer_OpenMP(res);

  767. 	if(data->adaptiveGrid==1){

  768. 		psidata = data->grid->x;

  769. 	}else{

  770. 		psidata = data->uniformGrid;

  771. 	}

  772. 

  773. 	innerMassFractionsData =  data->innerMassFractions;

  774. 

  775. 	/* Grid stencil:*/

  776. 

  777. 	/*-------|---------*---------|---------*---------|-------*/

  778. 	/*-------|---------*---------|---------*---------|-------*/

  779. 	/*-------|---------*---------|---------*---------|-------*/

  780. 	/*-------m-------mhalf-------j-------phalf-------p-------*/

  781. 	/*-------|---------*---------|---------*---------|-------*/

  782. 	/*-------|---------*---------|---------*---------|-------*/

  783. 	/*-------|<=======dxm=======>|<=======dxp=======>|-------*/

  784. 	/*-------|---------*<======dxav=======>*---------|-------*/

  785. 	/*-------|<================dxpm=================>|-------*/

  786. 

  787. 	/* Various variables defined for book-keeping and storing previously

  788. 	 * calculated values:

  789. 	 * rho		: densities at points  m, mhalf, j, p, and phalf.

  790. 	 * area		: the matric at points m, mhalf, j, p, and phalf.

  791. 	 * m 		: exponent that determines geometry;

  792. 	 * lambda	: thermal conductivities at mhalf and phalf.

  793. 	 * mdot		: mass flow rate at m, j, and p.

  794. 	 * X		: mole fractions at j and p.

  795. 	 * YV		: diffusion fluxes at mhalf and phalf.

  796. 	 * Tgrad	: temperature gradient at mhalf and phalf.

  797. 	 * Tav		: average temperature between two points.

  798. 	 * Pav		: average pressure between two points.

  799. 	 * Yav		: average mass fractions between two points.

  800. 	 * Xgradhalf	: mole fraction gradient at j.

  801. 	 * Cpb		: mass based bulk specific heat.

  802. 	 * tranTerm	: transient terms.

  803. 	 * advTerm	: advection terms.

  804. 	 * diffTerm	: diffusion terms.

  805. 	 * srcTerm	: source terms.

  806. 	 */

  807. 

  808. 	double rhomhalf,  rhom, lambdamhalf, YVmhalf[nsp], 

  809. 		 rho, 

  810. 		 rhophalf, lambdaphalf, YVphalf[nsp],

  811. 		 Cpb, Cvb,  Cp[nsp], Cpl, dHvl, rhol,  wdot[nsp],   enthalpy[nsp], energy[nsp],

  812. 		 tranTerm, diffTerm,    srcTerm, advTerm,

  813. 		 area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq;

  814. 	

  815. 	/*Aliases for difference coefficients:*/

  816. 	double cendfm, cendfc, cendfp;

  817. 	/*Aliases for various grid spacings:*/

  818. 	double dpsip, dpsiav, dpsipm, dpsim, dpsimm;

  819. 	dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;

  820. 	//double mass, mdotIn;

  821. 	double mass, massDrop;

  822. 	double sum, sum1, sum2, sum3;

  823. 

  824. 	size_t j,k;

  825. 	int m;

  826. 	m=data->metric;					//Unitless

  827. 	mass=data->mass;				//Units: kg

  828. 	//massDrop=data->massDrop;		//Units: kg

  829. 	//mdotIn=data->mdot*calc_area(R(npts),&m);	//Units: kg/s

  830. 

  831. //	/*evaluate properties at j=1*************************/

  832. 	setGas(data,ydata,1);

  833. 	rhom=data->gas->density();                                       

  834.     Cpb=data->gas->cp_mass();       	//J/kg/K

  835.     Cvb=data->gas->cv_mass();       	//J/kg/K

  836. 	//TODO: Create user input model for these. They should not be constant

  837. 	//Cpl=4.182e03;						//J/kg/K for liquid water

  838. 	//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  

  839.     //     - 9.62429538e+01*T(1)        + 1.27131046e+04;

  840.     

  841.     /*Update specific heat:Cpl for liquid n-heptane*/

  842.    /*Based on the existiong data,a linear expression is used*/

  843.    // Cpl= 12.10476*T(1) - 746.60143; // Unit:J/(kg*K), Need to be improved later

  844.     Cpl = 2242.871642; //Unit:J/(kg*K),range:25~520K,Ginnings&Furukawa,NIST

  845. 

  846. 	//dHvl=2.260e6;						//J/kg	heat of vaporization of water

  847. 	//double Tr = T(1)/647.1;				//Reduced Temperature: Droplet Temperature divided by critical temperature of water

  848. 	//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

  849. 	

  850.     double Tr = T(1)/540.2;            //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST

  851.     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

  852.     

  853.     /*Following section is related to the property of water*/

  854.     //rhol = 997.0;	

  855. 	//massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input

  856.     

  857.     /*Following section is related to the property of liquid n-heptane (mainly density rho)*/

  858.     rhol = 684.0;                          //Unit:kg/m^3

  859.     massDrop = 1.0/3.0*R(1)*R(1)*R(1)*rhol;  //Unit:kg, this is the mass of liquid droplet  

  860. 

  861.     aream= calc_area(R(1),&m);	    

  862.     

  863. 

  864. 	/*******************************************************************/

  865. 	/*Calculate values at j=2's m and mhalf*****************************/

  866. 

  867. 	getTransport(data, ydata, 1, &rhomhalf,&lambdamhalf,YVmhalf);

  868. 	areamhalf= calc_area(HALF*(R(1)+R(2)),&m);                                 

  869. 	aream = calc_area(R(1),&m);                                 

  870. 	areamhalfsq= areamhalf*areamhalf;                                   

  871.                                                                              

  872. 	/*******************************************************************/

  873. 

  874. 	/*Fill up res with left side (center) boundary conditions:**********/

  875. 	/*We impose zero fluxes at the center:*/                        

  876. 

  877. 

  878. 	/*Mass:*/                                                            

  879. 	if(data->quasiSteady){

  880. 		Rres(1)=Rdot(1); // Should remain satisfied for quasi-steady droplet

  881. 	}else{

  882. 		Rres(1)=Rdot(1) + Mdot(1)/rhol/aream; 

  883. 	}

  884.                                                                         

  885. 	/*Species:*/                                                    

  886. 	sum=ZERO;

  887. 	

  888. 	//TODO: Antoine Parameters should be part of user input, not hardcoded

  889. 	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e5*pow(10,4.6543-(1435.264/(T(1)-64.848)))*data->gas->molecularWeight(k_drop-1) 

  890. 	//			    / P(1) / data->gas->meanMolecularWeight();

  891. 	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)))*data->gas->molecularWeight(k_drop-1) 

  892. 	//			    / P(1) / data->gas->meanMolecularWeight();

  893. 	//Yres(1,k_drop)=Y(1,k_drop)-2.566785e-02;

  894. 	

  895.     /*Following using the Antoine Parameters to calculate the pressure of volatile component(n-heptane)*/

  896.     //Antoine parameter from NIST website

  897.     double p_i=0.0 ; 

  898.     p_i = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199))) ;    //unit:Pa,FOR N-HEPTANE 

  899.     //p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)));     //Unit:Pa,FOR WATER

  900.     Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) 

  901. 				    / P(1) / data->gas->meanMolecularWeight();

  902. 	

  903.     sum=sum+Y(1,k_drop);

  904. 

  905. 	Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass

  906. 

  907. 	for (k = 1; k <=nsp; k++) {                                        

  908. 		if(k!=k_bath && k!=k_drop){

  909. 			//TODO: May need to include chemical source term

  910. 			Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf;

  911. 			//Yres(1,k)=YVmhalf(k)*areamhalf;

  912. 		    sum=sum+Y(1,k);

  913. 			//if(fabs(Mdot(1))>1e-14){

  914. 		      		///Yres(1,k)=innerMassFractionsData[k-1]-

  915. 					///  Y(1,k)-

  916. 					///  (YVmhalf(k)*areamhalf)/Mdot(1);

  917. 					//Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf;

  918. 			//}

  919. 			//else{

  920. 		    //  		//Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];

  921. 		    //  		Yres(1,k)=Y(2,k)-Y(1,k);

  922. 			//	/*The below flux boundary condition makes the

  923. 			//	 * problem more prone to diverge. How does one

  924. 			//	 * fix this?*/

  925. 		    //  		//Yres(1,k)=YVmhalf(k);

  926. 			//}

  927. 		}

  928. 	}                                                                 

  929. 	Yres(1,k_bath)=ONE-sum-Y(1,k_bath);                                     

  930. 	

  931. 	/*Energy:*/                                                     

  932. 	if(data->dirichletInner){

  933. 		//Tres(1)=Tdot(1);

  934. 		//Tres(1)=T(1)-data->innerTemperature;

  935. 		Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHvl - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop));

  936. 	}else{

  937. 		//Tres(1)=Tdot(1) - 

  938. 		//	(rhomhalf*areamhalfsq*lambdamhalf*(T(2)-T(1))/((psi(2)-psi(1))*mass) - Mdot(1)*dHvl)

  939. 		//	/ (massDrop * Cpl);

  940. 		Tres(1)=Tdot(1) - 

  941. 			(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl)

  942. 			/ (massDrop * Cpl);

  943. 

  944. 

  945. 

  946. 		//Tres(1)=T(2)-T(1);

  947. 		//Tres(1)=Tdot(1) - (Pdot(1)/(rhom*Cpb))

  948. 		//	+(double)(data->metric+1)*(rhomhalf*lambdamhalf*areamhalfsq*(T(2)-T(1))/psi(2)-psi(1));

  949. 	}

  950. 

  951. 	/*Pressure:*/                                                        

  952. 	Pres(1)=P(2)-P(1);

  953. 	/*Fill up res with governing equations at inner points:*************/

  954. 	for (j = 2; j < npts; j++) {                                         

  955.                                                                              

  956. 		/*evaluate various mesh differences*///                      

  957.         	dpsip =        (psi(j+1) - psi(j)  )*mass;                            

  958.         	dpsim =        (psi(j)   - psi(j-1))*mass;                            

  959.         	dpsiav =  HALF*(psi(j+1) - psi(j-1))*mass;                            

  960.         	dpsipm =       (psi(j+1) - psi(j-1))*mass;                            

  961. 		/***********************************///                      

  962.                                                                              

  963. 		/*evaluate various central difference coefficients*/       

  964.         	cendfm = - dpsip / (dpsim*dpsipm);                               

  965.         	cendfc =   (dpsip-dpsim) / (dpsip*dpsim);                         

  966.         	cendfp =   dpsim / (dpsip*dpsipm);                               

  967. 		/**************************************************/

  968.                                                                    

  969.                                                                   

  970. 		/*evaluate properties at j*************************/

  971. 		setGas(data,ydata,j);

  972. 		rho=data->gas->density();		//kg/m^3

  973.         	Cpb=data->gas->cp_mass();       	//J/kg/K

  974.         	Cvb=data->gas->cv_mass();       	//J/kg/K

  975. 		data->gas->getNetProductionRates(wdot); //kmol/m^3

  976. 		data->gas->getEnthalpy_RT(enthalpy);	//unitless

  977. 		data->gas->getCp_R(Cp);			//unitless

  978. 		area = calc_area(R(j),&m);	        //m^2

  979.                                                                            

  980. 		/*evaluate properties at p*************************/

  981. 		getTransport(data, ydata, j, &rhophalf,&lambdaphalf,YVphalf);

  982. 		areaphalf= calc_area(HALF*(R(j)+R(j+1)),&m);	    

  983. 		areaphalfsq= areaphalf*areaphalf;		    

  984. 		/**************************************************///     

  985. 

  986. 		/*Evaporating Mass*/

  987. 		Mdotres(j)=Mdot(j) - Mdot(j-1); 

  988. 

  989. 		/*Mass:*/

  990. 		/* ∂r/∂ψ = 1/ρA */

  991. 		Rres(j)=((R(j)-R(j-1))/dpsim)-(TWO/(rhom*aream+rho*area));

  992. 

  993. 		/*Energy:*/

  994. 		/* ∂T/∂t = - ṁ(∂T/∂ψ) 

  995. 		 * 	   + (∂/∂ψ)(λρA²∂T/∂ψ) 

  996. 		 * 	   - (A/cₚ) ∑ YᵢVᵢcₚᵢ(∂T/∂ψ) 

  997. 		 * 	   - (1/ρcₚ)∑ ώᵢhᵢ 

  998. 		 * 	   + (1/ρcₚ)(∂P/∂t) */

  999. 		/*Notes: 

  1000. 		 * λ has units J/m/s/K.

  1001. 		 * YᵢVᵢ has units kg/m^2/s.

  1002. 		 * hᵢ has units J/kmol, so we must multiply the enthalpy

  1003. 		 * defined above (getEnthalpy_RT) by T (K) and the gas constant

  1004. 		 * (J/kmol/K) to get the right units.

  1005. 		 * cₚᵢ has units J/kg/K, so we must multiply the specific heat

  1006. 		 * defined above (getCp_R) by the gas constant (J/kmol/K) and

  1007. 		 * divide by the molecular weight (kg/kmol) to get the right

  1008. 		 * units.

  1009. 		 * */

  1010. 

  1011. 		//enthalpy formulation:

  1012. 		//tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));

  1013. 		tranTerm = Tdot(j);

  1014. 		sum=ZERO;

  1015. 		sum1=ZERO;

  1016. 		for (k = 1; k <=nsp; k++) {

  1017. 			sum=sum+wdot(k)*enthalpy(k);

  1018. 			sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho

  1019. 			     *HALF*(YVmhalf(k)+YVphalf(k));

  1020. 		}

  1021. 		sum=sum*Cantera::GasConstant*T(j);

  1022. 		sum1=sum1*Cantera::GasConstant;

  1023. 		diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)

  1024. 			      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )

  1025. 			     /(dpsiav*Cpb) )

  1026. 			  -(sum1*area*(cendfp*T(j+1)

  1027. 			              +cendfc*T(j)

  1028. 			              +cendfm*T(j-1))/Cpb);

  1029. 		srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cpb); // Qdot is forced heating to cause ignition, sum is the conversion of chemical enthalpy to sensible enthalpy

  1030. 		advTerm   = (Mdot(j)*(T(j)-T(j-1))/dpsim);

  1031. 		Tres(j)= tranTerm - (Pdot(j)/(rho*Cpb))

  1032. 			+advTerm 

  1033. 			-diffTerm

  1034. 			+srcTerm;

  1035. 

  1036. 	//	//energy formulation:

  1037. 	//	tranTerm = Tdot(j);

  1038. 	//	sum=ZERO;

  1039. 	//	sum1=ZERO;

  1040. 	//	sum2=ZERO;

  1041. 	//	sum3=ZERO;

  1042. 	//	for (k = 1; k <=nsp; k++) {

  1043. 	//		energy(k)=enthalpy(k)-ONE;

  1044. 	//		sum=sum+wdot(k)*energy(k);

  1045. 	//		sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho

  1046. 	//		     *HALF*(YVmhalf(k)+YVphalf(k));

  1047. 	//		sum2=sum2+(YVmhalf(k)/data->gas->molecularWeight(k-1));

  1048. 	//		sum3=sum3+(YVphalf(k)/data->gas->molecularWeight(k-1));

  1049. 	//	}

  1050. 	//	sum=sum*Cantera::GasConstant*T(j);

  1051. 	//	sum1=sum1*Cantera::GasConstant;

  1052. 	//	diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)

  1053. 	//		      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )

  1054. 	//		     /(dpsiav*Cvb) )

  1055. 	//		  -(sum1*area*(cendfp*T(j+1)

  1056. 	//		              +cendfc*T(j)

  1057. 	//		              +cendfm*T(j-1))/Cvb);

  1058. 	//	srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cvb);

  1059. 	//	advTerm   = (mdotIn*(T(j)-T(j-1))/dpsim);

  1060. 	//	advTerm = advTerm + (Cantera::GasConstant*T(j)*area/Cvb)*((sum3-sum2)/dpsiav);

  1061. 	//	Tres(j)= tranTerm

  1062. 	//		+advTerm 

  1063. 	//		-diffTerm

  1064. 	//		+srcTerm;

  1065. 

  1066. 		/*Species:*/

  1067. 		/* ∂Yᵢ/∂t = - ṁ(∂Yᵢ/∂ψ) 

  1068. 		 * 	    - (∂/∂ψ)(AYᵢVᵢ) 

  1069. 		 * 	    + (ώᵢWᵢ/ρ)  */

  1070. 

  1071. 		sum=ZERO;

  1072. 		for (k = 1; k <=nsp; k++) {

  1073. 			if(k!=k_bath){

  1074. 				tranTerm  = Ydot(j,k);

  1075. 				diffTerm  = (YVphalf(k)*areaphalf

  1076. 					    -YVmhalf(k)*areamhalf)/dpsiav;

  1077. 				srcTerm   = wdot(k)

  1078. 			     		    *(data->gas->molecularWeight(k-1))/rho;

  1079. 				advTerm   = (Mdot(j)*(Y(j,k)-Y(j-1,k))/dpsim);

  1080. 

  1081. 				Yres(j,k)= tranTerm 

  1082. 					  +advTerm

  1083. 					  +diffTerm

  1084. 					  -srcTerm;

  1085. 

  1086. 		      		sum=sum+Y(j,k);

  1087. 			}

  1088. 		}

  1089. 		Yres(j,k_bath)=ONE-sum-Y(j,k_bath);

  1090. 

  1091. 		/*Pressure:*/

  1092. 		Pres(j) = P(j+1)-P(j);

  1093. 

  1094. 		/*Assign values evaluated at p and phalf to m       

  1095. 		 * and mhalf to save some cpu cost:****************/

  1096. 		areamhalf=areaphalf;                                

  1097. 		areamhalfsq=areaphalfsq;                            

  1098. 		aream=area;

  1099. 		rhom=rho;

  1100. 		rhomhalf=rhophalf;                                  

  1101. 		lambdamhalf=lambdaphalf;                            

  1102. 		for (k = 1; k <=nsp; k++) {                         

  1103. 			YVmhalf(k)=YVphalf(k);                      

  1104. 		}                                                   

  1105. 		/**************************************************/

  1106. 		                                                    

  1107. 	}                                                           

  1108. 	/*******************************************************************///

  1109. 

  1110. 

  1111. 	/*Fill up res with right side (wall) boundary conditions:***********/

  1112. 	/*We impose zero fluxes at the wall:*/                               

  1113. 	setGas(data,ydata,npts);

  1114. 	rho=data->gas->density();                               

  1115. 	area = calc_area(R(npts),&m);				   

  1116. 

  1117. 	/*Mass:*/                                                          

  1118. 	dpsim=(psi(npts)-psi(npts-1))*mass;

  1119. 	Rres(npts)=((R(npts)-R(npts-1))/dpsim)-(TWO/(rhom*aream+rho*area));

  1120. 

  1121. 	/*Energy:*/                                             

  1122. 	if(data->dirichletOuter){

  1123. 		//Tres(npts)=T(npts)-data->wallTemperature;

  1124. 		Tres(npts)=Tdot(npts);

  1125. 	}

  1126. 	else{

  1127. 		Tres(npts)=T(npts)-T(npts-1);

  1128. 	}

  1129.                                                                 

  1130. 	/*Species:*/                                            

  1131. 	sum=ZERO;

  1132. 	if(data->dirichletOuter){

  1133. 		for (k = 1; k <=nsp; k++) {                             

  1134. 			if(k!=k_bath){

  1135. 			    	Yres(npts,k)=Ydot(npts,k);

  1136. 			    	sum=sum+Y(npts,k);

  1137. 			}

  1138. 		}                                                       

  1139. 	}

  1140. 	else{

  1141. 		for (k = 1; k <=nsp; k++) {                             

  1142. 			if(k!=k_bath){

  1143. 			    	Yres(npts,k)=Y(npts,k)-Y(npts-1,k);

  1144. 			    	//Yres(npts,k)=YVmhalf(k);         

  1145. 			    	sum=sum+Y(npts,k);

  1146. 			}

  1147. 		}                                                       

  1148. 	}

  1149. 	Yres(npts,k_bath)=ONE-sum-Y(npts,k_bath);             

  1150.                                                                 

  1151. 	

  1152. 	/*Pressure:*/                                                      

  1153. 	if(data->constantPressure){

  1154. 		Pres(npts)=Pdot(npts)-data->dPdt;                                     

  1155. 	}

  1156. 	else{

  1157. 		Pres(npts)=R(npts)-data->domainLength;                     

  1158. 		//Pres(npts)=Rdot(npts);                     

  1159. 	}

  1160. 

  1161. 	/*Evaporating Mass*/

  1162. 	Mdotres(npts)=Mdot(npts) - Mdot(npts-1); 

  1163. 

  1164. 	//for (j = 1; j <=npts; j++) {

  1165. 	//	//for (k = 1; k <=nsp; k++) {

  1166. 	//	//	Yres(j,k)=Ydot(j,k);

  1167. 	//	//}

  1168. 	//	//Tres(j)=Tdot(j);

  1169. 	//}

  1170. 

  1171. 	return(0);

  1172. 

  1173. }

  1174. 

  1175. void printSpaceTimeHeader(UserData data)

  1176. {

  1177. 	fprintf((data->output), "%15s\t","#1");

  1178. 	for (size_t k = 1; k <=data->nvar+1; k++) {

  1179. 		fprintf((data->output), "%15lu\t",k+1);

  1180. 	}

  1181. 	fprintf((data->output), "%15lu\n",data->nvar+3);

  1182. 

  1183. 	fprintf((data->output), "%15s\t%15s\t%15s\t","#psi","time(s)","dpsi");

  1184. 	fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");

  1185. 	for (size_t k = 1; k <=data->nsp; k++) {

  1186. 		fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());

  1187. 	}

  1188. 	fprintf((data->output), "%15s\t","Pressure(Pa)");

  1189. 	fprintf((data->output), "%15s\n","Mdot (kg/s)");

  1190. }

  1191. 

  1192. void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)

  1193. {

  1194. 	double *ydata,*psidata;

  1195.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1196. 

  1197. 	if(data->adaptiveGrid){

  1198.   		psidata  = data->grid->x;

  1199. 	}else{

  1200.   		psidata  = data->uniformGrid;

  1201. 	}

  1202. 

  1203. 	for (size_t i = 0; i < data->npts; i++) {

  1204. 		fprintf(output, "%15.9e\t%15.9e\t",psi(i+1),t);

  1205. 		if(i==0){

  1206. 			fprintf(output, "%15.9e\t",psi(2)-psi(1));

  1207. 		}

  1208. 		else{

  1209. 			fprintf(output, "%15.6e\t",psi(i+1)-psi(i));

  1210. 		}

  1211. 		for (size_t j = 0; j < data->nvar; j++) {

  1212. 			if(j!=data->nvar-1){

  1213. 				fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);

  1214. 			}else{

  1215. 				fprintf(output, "%15.9e",ydata[j+i*data->nvar]);

  1216. 			}

  1217. 		}

  1218. 		fprintf(output, "\n");

  1219. 	}

  1220. 	fprintf(output, "\n");

  1221. }

  1222. 

  1223. void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data){

  1224. 	double *ydata,*psidata, *ydotdata;

  1225. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1226. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1227. 	if(data->adaptiveGrid){

  1228.   		psidata  = data->grid->x;

  1229. 	}else{

  1230.   		psidata  = data->uniformGrid;

  1231. 	}

  1232.  	fwrite(&t, sizeof(t), 1, output);		//write time

  1233.  	fwrite(psidata, data->npts*sizeof(psidata), 1, output);	//write grid

  1234.  	fwrite(ydata, data->neq*sizeof(ydata), 1, output);	//write solution

  1235.  	fwrite(ydotdata, data->neq*sizeof(ydotdata), 1, output);	//write solutiondot

  1236. }

  1237. 

  1238. void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data){

  1239. 	double *ydata,*psidata, *ydotdata;

  1240. 	double t;

  1241. 	if(data->adaptiveGrid){

  1242.   		psidata  = data->grid->x;

  1243. 	}else{

  1244.   		psidata  = data->uniformGrid;

  1245. 	}

  1246. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1247. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1248.  	fread(&t, sizeof(t), 1, input);

  1249. 	data->tNow=t;

  1250.  	fread(psidata, data->npts*sizeof(psidata), 1, input);

  1251.  	fread(ydata, data->neq*sizeof(ydata), 1, input);

  1252.  	fread(ydotdata, data->neq*sizeof(ydotdata), 1, input);

  1253. 	if(data->adaptiveGrid){

  1254. 		storeGrid(data->grid->x,data->grid->xOld,data->npts);

  1255. 	}

  1256. }

  1257. 

  1258. void printGlobalHeader(UserData data)

  1259. {

  1260. 	fprintf((data->globalOutput), "%8s\t","#Time(s)");

  1261. 	//fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");

  1262. 	fprintf((data->globalOutput), "%15s","bFlux(kg/m^2/s)");

  1263. 	fprintf((data->globalOutput), "%16s","  IsothermPos(m)");

  1264. 	fprintf((data->globalOutput), "%15s\t","Pressure(Pa)");

  1265. 	fprintf((data->globalOutput), "%15s\t","Pdot(Pa/s)");

  1266. 	fprintf((data->globalOutput), "%15s\t","gamma");

  1267. 	fprintf((data->globalOutput), "%15s\t","S_u(m/s)");

  1268. 	fprintf((data->globalOutput), "%15s\t","Tu(K)");

  1269. 	fprintf((data->globalOutput), "\n");

  1270. }

  1271. //

  1272. //void printSpaceTimeRates(double t, N_Vector ydot, UserData data)

  1273. //{

  1274. //	double *ydotdata,*psidata;

  1275. //  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);

  1276. //  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);

  1277. //	for (int i = 0; i < data->npts; i++) {

  1278. //		fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);

  1279. //		for (int j = 0; j < data->nvar; j++) {

  1280. //			fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);

  1281. //		}

  1282. //		fprintf((data->ratesOutput), "\n");

  1283. //	}

  1284. //	fprintf((data->ratesOutput), "\n\n");

  1285. //}

  1286. //

  1287. void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)

  1288. {

  1289. 	double *ydata,*ydotdata, *innerMassFractionsData, *psidata;

  1290. 	innerMassFractionsData =  data->innerMassFractions;

  1291. 

  1292. 	if(data->adaptiveGrid){

  1293.   		psidata  = data->grid->x;

  1294. 	}else{

  1295.   		psidata  = data->uniformGrid;

  1296. 	}

  1297. 

  1298. 	double TAvg, RAvg, YAvg, psiAvg;

  1299.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1300.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1301. 	TAvg=data->isotherm;

  1302. 	double sum=ZERO;

  1303. 	double dpsim,area,aream,drdt;

  1304. 	double Cpb,Cvb,gamma,rho,flameArea,Tu;

  1305. 

  1306. 

  1307. 	/*Find the isotherm chosen by the user*/

  1308. 	size_t j=1;

  1309. 	size_t jj=1;

  1310. 	size_t jjj=1;

  1311. 	double wdot[data->nsp];

  1312. 	double wdotMax=0.0e0;

  1313. 	double advTerm=0.0e0;

  1314. 	psiAvg=0.0e0;

  1315. 	if(T(data->npts)>T(1)){

  1316. 		while (T(j)<TAvg) {

  1317. 			j=j+1;

  1318. 		}

  1319. 		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(data->npts,data->k_oxidizer);

  1320. 		while (fabs((T(jj+1)-T(jj))/T(jj))>1e-08) {

  1321. 			jj=jj+1;

  1322. 		}

  1323. 

  1324. 		setGas(data,ydata,jj);

  1325. 		Tu=T(jj);

  1326. 		rho=data->gas->density();

  1327.         	Cpb=data->gas->cp_mass();       	//J/kg/K

  1328.         	Cvb=data->gas->cv_mass();       	//J/kg/K

  1329. 		gamma=Cpb/Cvb;

  1330. 		

  1331. 	}

  1332. 	else{

  1333. 		while (T(j)>TAvg) {

  1334. 			j=j+1;

  1335. 		}

  1336. 		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(1,data->k_oxidizer);

  1337. 		while (fabs((T(data->npts-jj-1)-T(data->npts-jj))/T(data->npts-jj))>1e-08) {

  1338. 			jj=jj+1;

  1339. 		}

  1340. 

  1341. 		setGas(data,ydata,data->npts-jj);

  1342. 		Tu=T(data->npts-jj);

  1343. 		rho=data->gas->density();

  1344.         	Cpb=data->gas->cp_mass();       	//J/kg/K

  1345.         	Cvb=data->gas->cv_mass();       	//J/kg/K

  1346. 		gamma=Cpb/Cvb;

  1347. 	}

  1348. 	if(T(j)<TAvg){

  1349. 		RAvg=((R(j+1)-R(j))/(T(j+1)-T(j)))*(TAvg-T(j))+R(j);

  1350. 	}

  1351. 	else{

  1352. 		RAvg=((R(j)-R(j-1))/(T(j)-T(j-1)))*(TAvg-T(j-1))+R(j-1);

  1353. 	}

  1354. 

  1355. 	////Experimental burning speed calculation:	

  1356. 	//int nMax=0;

  1357. 	////nMax=maxCurvIndex(ydata, data->nt, data->nvar,

  1358. 	////     	data->grid->x, data->npts);

  1359. 	//nMax=maxGradIndex(ydata, data->nt, data->nvar,

  1360. 	//     	data->grid->x, data->npts);

  1361. 	//advTerm=(T(nMax)-T(nMax-1))/(data->mass*(psi(nMax)-psi(nMax-1)));

  1362. 	//aream=calc_area(R(nMax),&data->metric);	    

  1363. 	////setGas(data,ydata,nMax);

  1364. 	////rho=data->gas->density();

  1365. 	//psiAvg=-Tdot(nMax)/(rho*aream*advTerm);

  1366. 	////if(t>data->ignTime){

  1367. 	////	for(size_t n=2;n<data->npts;n++){

  1368. 	////		setGas(data,ydata,n);

  1369. 	////		data->gas->getNetProductionRates(wdot); //kmol/m^3

  1370. 	////		advTerm=(T(n)-T(n-1))/(data->mass*(psi(n)-psi(n-1)));

  1371. 	////		if(fabs(wdot[data->k_oxidizer-1])>=wdotMax){

  1372. 	////			aream=calc_area(R(n),&data->metric);	    

  1373. 	////			psiAvg=-Tdot(n)/(rho*aream*advTerm);

  1374. 	////			wdotMax=fabs(wdot[data->k_oxidizer-1]);

  1375. 	////		}

  1376. 	////	}

  1377. 	////}

  1378. 	////else{

  1379. 	////	psiAvg=0.0e0;

  1380. 	////}

  1381. 

  1382. 

  1383. 	//drdt=(RAvg-data->flamePosition[1])/(t-data->flameTime[1]);

  1384. 	//data->flamePosition[0]=data->flamePosition[1];

  1385. 	//data->flamePosition[1]=RAvg;

  1386. 	//data->flameTime[0]=data->flameTime[1];

  1387. 	//data->flameTime[1]=t;

  1388. 	//flameArea=calc_area(RAvg,&data->metric);

  1389. 

  1390. 	/*Use the Trapezoidal rule to calculate the mass burning rate based on

  1391. 	 * the consumption of O2*/

  1392. 	aream= calc_area(R(1)+1e-03*data->domainLength,&data->metric);	    

  1393. 	for (j = 2; j <data->npts; j++) {

  1394.         	dpsim=(psi(j)-psi(j-1))*data->mass;                            

  1395. 		area= calc_area(R(j),&data->metric);	    

  1396. 		sum=sum+HALF*dpsim*((Ydot(j-1,data->k_oxidizer)/aream)

  1397. 				   +(Ydot(j,data->k_oxidizer)/area));

  1398. 		aream=area;

  1399. 	}

  1400. 

  1401. 	//double maxOH,maxHO2;

  1402. 	//maxOH=0.0e0;

  1403. 	//maxHO2=0.0e0;

  1404. 	//for(j=1;j<data->npts;j++){

  1405. 	//	if(Y(j,data->k_OH)>maxOH){

  1406. 	//		maxOH=Y(j,data->k_OH);

  1407. 	//	}

  1408. 	//}

  1409. 	//for(j=1;j<data->npts;j++){

  1410. 	//	if(Y(j,data->k_HO2)>maxHO2){

  1411. 	//		maxHO2=Y(j,data->k_HO2);

  1412. 	//	}

  1413. 	//}

  1414. 

  1415. 	fprintf((data->globalOutput), "%15.6e\t",t);

  1416. 	//fprintf((data->globalOutput), "%15.6e\t",psiAvg);

  1417. 	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/YAvg);

  1418. 	fprintf((data->globalOutput), "%15.6e\t",RAvg);

  1419. 	fprintf((data->globalOutput), "%15.6e\t",P(data->npts));

  1420. 	fprintf((data->globalOutput), "%15.6e\t",Pdot(data->npts));

  1421. 	fprintf((data->globalOutput), "%15.6e\t",gamma);

  1422. 	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/(YAvg*rho));

  1423. 	fprintf((data->globalOutput), "%15.6e\t",Tu);

  1424. 	fprintf((data->globalOutput), "\n");

  1425. }

  1426. 

  1427. //

  1428. //void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data)

  1429. //{

  1430. //	double *ydata,*psidata;

  1431. //  	ydata    = N_VGetArrayPointer_OpenMP(y);

  1432. //  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);

  1433. //	for (int i = 0; i < data->npts; i++) {

  1434. //		fprintf((data->gridOutput), "%15.6e\t%15.6e\t",psi(i+1),t);

  1435. //		for (int j = 0; j < data->nvar; j++) {

  1436. //			fprintf((data->gridOutput), "%15.6e\t",ydata[j+i*data->nvar]);

  1437. //		}

  1438. //		fprintf((data->gridOutput), "\n");

  1439. //	}

  1440. //	fprintf((data->gridOutput), "\n\n");

  1441. //}

  1442. //

  1443. //

  1444. ////void repairSolution(N_Vector y, N_Vector ydot, UserData data){

  1445. ////	int npts=data->npts;

  1446. ////  	double *ydata;

  1447. ////  	double *ydotdata;

  1448. ////	ydata    = N_VGetArrayPointer_OpenMP(y); 

  1449. ////	ydotdata    = N_VGetArrayPointer_OpenMP(ydot); 

  1450. ////

  1451. ////	T(2)=T(1);

  1452. ////	T(npts-1)=T(npts);

  1453. ////	Tdot(2)=Tdot(1);

  1454. ////	Tdot(npts-1)=Tdot(npts);

  1455. ////	for (int k = 1; k <=data->nsp; k++) {

  1456. ////		    Y(2,k)=Y(1,k);

  1457. ////		    Y(npts-1,k)=Y(npts,k);

  1458. ////

  1459. ////		    Ydot(2,k)=Ydot(1,k);

  1460. ////		    Ydot(npts-1,k)=Ydot(npts,k);

  1461. ////	}

  1462. ////}

  1463. 

  1464. 

  1465. /*Following functions are added to derive the characteristic time scale of species*/

  1466. void getTimescale(UserData data, N_Vector* y){

  1467. 	size_t i, k, nsp,npts ; 

  1468. 	nsp = data->nsp ; 

  1469. 	npts = data->npts ; 

  1470. 	double rho, wdot_mole[nsp], wdot_mass[nsp],MW[nsp],concentra[nsp]; 

  1471. 	//double time_scale[npts*nsp] ;   

  1472. 

  1473. 	double* ydata ; 

  1474. 	ydata = N_VGetArrayPointer_OpenMP(*y); //access the data stored in N_Vector* y 

  1475. 	

  1476. 	for(i=1; i<= npts;i++){

  1477. 		setGas(data,ydata,i); //set the gas state at each grid point

  1478. 		rho = data->gas->density(); //get the averaged density at each grid point, Unit:kg/m^3

  1479. 		data->gas->getNetProductionRates(wdot_mole) ; //Unit:kmol/(m^3 * s)

  1480. 		

  1481. 		for(k=1; k<= nsp; k++){

  1482. 			MW(k) = data->gas->molecularWeight(k-1) ; //Unit:kg/kmol

  1483. 		}

  1484. 

  1485. 		for(k=1; k<= nsp; k++){

  1486. 			wdot_mass(k) = wdot_mole(k) * MW(k) ; //Unit:kg/(m^3 * s)

  1487. 		}

  1488. 

  1489. 		for(k=1; k<= nsp; k++){

  1490. 			concentra(k) = Y(i,k) * rho ;  //Unit:kg/m^3

  1491. 		}

  1492. 

  1493. 		for(k=1;k<= nsp;k++){

  1494. 			data->time_scale(i,k) = concentra(k)/(wdot_mass(k)+ 1.00e-16) ; 

  1495. 		}

  1496.         

  1497. 

  1498. 

  1499. 	}

  1500.     

  1501. }

  1502. 

  1503. 

  1504. void printTimescaleHeader(UserData data)

  1505. {

  1506. 	fprintf((data->timescaleOutput), "%15s\t","#1");

  1507. 	for (size_t k = 1; k <=data->nsp+1; k++) {

  1508. 		fprintf((data->timescaleOutput), "%15lu\t",k+1);

  1509. 	}

  1510. 	fprintf((data->timescaleOutput), "%15lu\n",data->nsp+3);

  1511. 

  1512. 

  1513. 	fprintf((data->timescaleOutput), "%15s\t%15s\t%15s\t","#time","radius","Temp(K)");

  1514. 	//fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");

  1515. 	for (size_t k = 1; k <=(data->nsp); k++) {

  1516. 		fprintf((data->timescaleOutput), "%15s\t",data->gas->speciesName(k-1).c_str());

  1517. 	}

  1518. 	//fprintf((data->output), "%15s\t","Pressure(Pa)");

  1519. 	//fprintf((data->timescaleOutput), "%15s\n",data->gas->speciesName(data->nsp-1).c_str());

  1520. 	//fprintf((data->output), "%15s\n","Mdot (kg/s)");

  1521. 	fprintf((data->timescaleOutput), "\n");

  1522. }

  1523. 

  1524. void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data)

  1525. {

  1526. 	double* ydata; 

  1527. 	ydata = N_VGetArrayPointer_OpenMP(*y); 

  1528. 

  1529. 	for(size_t i=1 ; i<= data->npts; i++){

  1530. 		fprintf(output, "%15.9e\t%15.9e\t",t,R(i));

  1531. 		fprintf(output, "%15.9e\t",T(i));

  1532.         

  1533. 		for(size_t k=1; k<=data->nsp;k++){

  1534. 			fprintf(output, "%15.9e\t",data->time_scale(i,k));

  1535. 		}

  1536. 		fprintf(output, "\n");

  1537. 	}

  1538. 	fprintf(output, "\n");

  1539. 

  1540. }

  1541. 

  1542. void floorSmallValue(UserData data, N_Vector* y)

  1543. {

  1544.     double* ydata;

  1545.     ydata = N_VGetArrayPointer_OpenMP(*y);

  1546.     //double sum = 0.00; 

  1547.     size_t k_bath = data->k_bath;

  1548.    

  1549. 	/*Floor small values to zero*/

  1550. 	for (size_t i = 1; i <=data->npts; i++) {

  1551. 		for (size_t k = 1; k <=data->nsp; k++) {

  1552. 			if(fabs(Y(i,k))<=data->massFractionTolerance){

  1553. 				Y(i,k)=0.0e0;

  1554. 			}

  1555. 		}

  1556. 	}

  1557.     

  1558.     /*Dump the error to the bath gas*/

  1559. 	for (size_t i = 1; i <=data->npts; i++) {

  1560.         double sum = 0.00;

  1561. 		for (size_t k = 1; k <=data->nsp; k++) {

  1562. 			if(k!=k_bath){

  1563. 				sum = sum + Y(i,k);

  1564. 			}

  1565. 		}

  1566.         Y(i,k_bath) = ONE-sum;

  1567. 	}

  1568.     

  1569. 

  1570. }

  1571. 

  1572. void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv)

  1573. {

  1574.     double* ydata;

  1575.     realtype* atolvdata; 

  1576.     ydata = N_VGetArrayPointer_OpenMP(*y);

  1577.     atolvdata = N_VGetArrayPointer_OpenMP(*atolv);

  1578.     double relTol,radTol,tempTol,presTol,massFracTol,bathGasTol,mdotTol;

  1579.     double maxT=0.00, ignT=1800;

  1580.     relTol = 1.0e-03 ; 

  1581.     radTol = 1.0e-05 ; 

  1582.     tempTol = 1.0e-01 ; 

  1583.     presTol = 1.0e00 ; 

  1584.     massFracTol = 1.0e-08; 

  1585.     bathGasTol = 1.0e-04 ; 

  1586.     mdotTol = 1.0e-09 ; 

  1587. 

  1588. 

  1589.     /*Get the maximum Tempture*/

  1590.     for (size_t i =1;i <= data->npts;i++){

  1591.         if(T(i) > maxT){

  1592.             maxT = T(i);

  1593.         }

  1594.     }

  1595. 

  1596.     /*reset the tolerance when maxT > ignT*/

  1597.     if(maxT >= ignT){

  1598.         data->relativeTolerance = relTol;

  1599. 

  1600.         for(size_t i =1; i<=data->npts; i++){

  1601.             atolT(i) = tempTol; 

  1602.             atolR(i) = radTol ; 

  1603.             atolP(i) = presTol ;

  1604.             atolMdot(i) = mdotTol ; 

  1605.             

  1606.             for(size_t k =1; k<= data->nsp; k++){

  1607.                 if(k!=data->k_bath){

  1608.                     atolY(i,k) = massFracTol;

  1609.                 }else{

  1610.                     atolY(i,k) = bathGasTol;

  1611.                 }

  1612.             }

  1613.         }

  1614.     }

  1615. 

  1616. }