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.
7 Commits
1 Branch
1.7MB
Tree: 293673d25d
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from '293673d25d'
${ noResults }
Droplet_LTORC/residue.cpp

1762 lines
52KB
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 maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double *x ,size_t nPts){

  13.     double maxT = 0.0e0;

  14.     double TempT = 0.0e0;

  15.     double pos = 0.0e0;

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

  17.         TempT = y[(i-1)*nvar+nt] ; 

  18.         if(TempT > maxT){

  19.             maxT = TempT ; 

  20.             pos = x[i-1];

  21.         }

  22.     }

  23.     return(pos); 

  24. }

  25. 

  26. double maxTemperature(const double* y,const size_t nt,const size_t nvar ,size_t nPts){

  27.     double maxT = 0.0e0;

  28.     double TempT = 0.0e0;

  29.     //int index = 0 ; 

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

  31.         TempT = y[(i-1)*nvar+nt] ; 

  32.         if(TempT > maxT){

  33.             maxT = TempT ; 

  34.         }

  35.     }

  36.     return(maxT); 

  37. }

  38. 

  39. int maxTemperatureIndex(const double* y,const size_t nt,const size_t nvar ,size_t nPts){

  40.     double maxT = 0.0e0;

  41.     double TempT = 0.0e0;

  42.     int index = 0 ; 

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

  44.         TempT = y[(i-1)*nvar+nt] ; 

  45.         if(TempT > maxT){

  46.             maxT = TempT ; 

  47.             index = i; 

  48.         }

  49.     }

  50.     return(index); 

  51. }

  52. 

  53. double maxCurvPositionR(const double* y, const size_t nt, 

  54. 	               const size_t nvar, const size_t nr, size_t nPts){

  55. 	double maxCurvT=0.0e0;

  56. 	double gradTp=0.0e0;

  57. 	double gradTm=0.0e0;

  58. 	double curvT=0.0e0;

  59. 	double dx=0.0e0;

  60.     double dr=0.0e0;

  61. 	double pos=0.0e0;

  62. 	size_t j,jm,jp;

  63.     size_t r,rp,rm;

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

  65. 		j=i*nvar+nt;

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

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

  68.         r = i*nvar+nr;

  69.         rm = (i-1)*nvar+nr;

  70.         rp = (i+1)*nvar+nr; 

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

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

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

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

  75.         gradTp = fabs((y[jp]-y[j])/(y[rp]-y[r]));

  76.         gradTp = fabs((y[j]-y[jm])/(y[r]-y[rm]));

  77. 		dr = 0.5e0*(y[rp]-y[rm]);

  78.         curvT=(gradTp-gradTm)/dr; 

  79.         if (curvT>=maxCurvT) {

  80. 			maxCurvT=curvT;

  81. 			pos=y[r];

  82. 		}

  83. 	}

  84. 	return(pos);

  85. }

  86. 

  87. int maxCurvIndexR(const double* y, const size_t nt, 

  88. 	               const size_t nvar, const size_t nr, size_t nPts){

  89. 	double maxCurvT=0.0e0;

  90. 	double gradTp=0.0e0;

  91. 	double gradTm=0.0e0;

  92. 	double curvT=0.0e0;

  93. 	double dx=0.0e0;

  94.     double dr=0.0e0;

  95. 	int pos=0;

  96. 	size_t j,jm,jp;

  97.     size_t r,rm,rp;

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

  99. 		j=i*nvar+nt;

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

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

  102.         r = i*nvar+nr;

  103.         rm = (i-1)*nvar+nr;

  104.         rp = (i+1)*nvar+nr; 

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

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

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

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

  109.         gradTp = fabs((y[jp]-y[j])/(y[rp]-y[r]));

  110.         gradTp = fabs((y[j]-y[jm])/(y[r]-y[rm]));

  111. 		dr = 0.5e0*(y[rp]-y[rm]);

  112.         curvT=(gradTp-gradTm)/dr; 

  113.         if (curvT>=maxCurvT) {

  114. 			maxCurvT=curvT;

  115. 			pos=i;

  116. 		}

  117. 	}

  118. 	return(pos);

  119. }

  120. 

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

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

  123. 	double maxGradT=0.0e0;

  124. 	double gradT=0.0e0;

  125. 	double pos=0.0e0;

  126. 	size_t j,jm;

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

  128. 		j=i*nvar+nt;

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

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

  131. 		if (gradT>=maxGradT) {

  132. 			maxGradT=gradT;

  133. 			pos=x[i];

  134. 		}

  135. 	}

  136. 	return(pos);

  137. }

  138. 

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

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

  141. 	double maxGradT=0.0e0;

  142. 	double gradT=0.0e0;

  143. 	int pos=0;

  144. 	size_t j,jm;

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

  146. 		j=i*nvar+nt;

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

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

  149. 		if (gradT>=maxGradT) {

  150. 			maxGradT=gradT;

  151. 			pos=i;

  152. 		}

  153. 	}

  154. 	return(pos);

  155. }

  156. 

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

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

  159. 	double maxCurvT=0.0e0;

  160. 	double gradTp=0.0e0;

  161. 	double gradTm=0.0e0;

  162. 	double curvT=0.0e0;

  163. 	double dx=0.0e0;

  164. 	double pos=0.0e0;

  165. 	size_t j,jm,jp;

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

  167. 		j=i*nvar+nt;

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

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

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

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

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

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

  174. 		if (curvT>=maxCurvT) {

  175. 			maxCurvT=curvT;

  176. 			pos=x[i];

  177. 		}

  178. 	}

  179. 	return(pos);

  180. }

  181. 

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

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

  184. 	double maxCurvT=0.0e0;

  185. 	double gradTp=0.0e0;

  186. 	double gradTm=0.0e0;

  187. 	double curvT=0.0e0;

  188. 	double dx=0.0e0;

  189. 	int pos=0;

  190. 	size_t j,jm,jp;

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

  192. 		j=i*nvar+nt;

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

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

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

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

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

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

  199. 		if (curvT>=maxCurvT) {

  200. 			maxCurvT=curvT;

  201. 			pos=i;

  202. 		}

  203. 	}

  204. 	return(pos);

  205. }

  206. 

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

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

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

  210. //	size_t j;

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

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

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

  214. //			pos=x[i];

  215. //			break;

  216. //		}

  217. //	}

  218. //	return(pos);

  219. //}

  220. 

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

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

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

  224. 	double pos=x[0];

  225. 	size_t j;

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

  227. 		j=i*nvar+nt;

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

  229. 			pos=x[i];

  230. 			break;

  231. 		}

  232. 	}

  233. 	return(pos);

  234. }

  235. 

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

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

  238. 

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

  240. 

  241. 	gsl_interp_accel* acc;

  242. 	gsl_spline* spline;

  243. 	acc = gsl_interp_accel_alloc();

  244. 	spline = gsl_spline_alloc(gsl_interp_steffen, nPts);

  245. 

  246. 	gsl_interp_accel* accdot;

  247. 	gsl_spline* splinedot;

  248. 	accdot    = gsl_interp_accel_alloc();

  249. 	splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);

  250. 

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

  252. 

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

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

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

  256. 		}

  257. 

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

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

  260. 

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

  262. 

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

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

  265. 		}

  266. 	}

  267. 

  268. 	//Exploring "fixing" boundary conditions:

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

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

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

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

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

  274. 	//}

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

  276. 	

  277. 	gsl_interp_accel_free(acc);

  278. 	gsl_spline_free(spline);

  279. 

  280. 	gsl_interp_accel_free(accdot);

  281. 	gsl_spline_free(splinedot);

  282. }

  283. 

  284. //Locate bath gas:

  285. size_t BathGasIndex(UserData data){

  286. 	size_t index=0;

  287. 	double max1;

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

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

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

  291. 		if(max1>=max){

  292. 			max=max1;

  293. 			index=k;

  294. 		}

  295. 	}

  296. 	return(index+1);

  297. }

  298. 

  299. //Locate Oxidizer:

  300. size_t oxidizerIndex(UserData data){

  301. 	size_t index=0;

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

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

  304. 			index=k;

  305. 		}

  306. 	}

  307. 	return(index);

  308. }

  309. 

  310. //Locate OH:

  311. size_t OHIndex(UserData data){

  312. 	size_t index=0;

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

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

  315. 			index=k;

  316. 		}

  317. 	}

  318. 	return(index);

  319. }

  320. 

  321. //Locate HO2:

  322. size_t HO2Index(UserData data){

  323. 	size_t index=0;

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

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

  326. 			index=k;

  327. 		}

  328. 	}

  329. 	return(index);

  330. }

  331. 

  332. //Locate species index:

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

  334. 	size_t index=0;

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

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

  337. 			index=k;

  338. 		}

  339. 	}

  340. 	return(index);

  341. }

  342. 

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

  344. 	double *iddata;

  345.   	N_VConst(ONE, *id);

  346.   	iddata   = N_VGetArrayPointer_OpenMP(*id);

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

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

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

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

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

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

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

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

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

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

  357. 		Rid(i)=ZERO;

  358. 		Pid(i)=ZERO;

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

  360. 		Mdotid(i)=ZERO;

  361. 	}	

  362. 	Mdotid(1)=ZERO;

  363. 	Rid(1)=ONE;

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

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

  366. 	if(data->constantPressure){

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

  368. 	}else{

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

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

  371. 	}

  372. 	if(data->dirichletInner){

  373. 		Tid(1)=ZERO; 

  374. 	}else{

  375. 		Tid(1)=ONE;

  376. 	}

  377. 	if(data->dirichletOuter){

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

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

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

  381. 			}

  382. 			Yid(1,k)=ZERO;

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

  384. 		}

  385. 	}else{

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

  387. 			Yid(1,k)=ZERO;

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

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

  390. 		}

  391. 	}

  392. 	return(0);

  393. }

  394. 

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

  396. 	switch (*i) {

  397. 		case 0:

  398. 			return(ONE);

  399. 		case 1:

  400. 			return(x);

  401. 		case 2:

  402. 			return(x*x);

  403. 		default:

  404. 			return(ONE);

  405. 	}

  406. }

  407. 

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

  409. 

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

  411. 

  412. 	size_t bufLen=10000;

  413. 	size_t nRows=0;

  414. 	size_t nColumns=nvar;

  415. 

  416. 	char buf[bufLen];

  417. 	char buf1[bufLen];

  418. 	char comment[1];

  419. 	char *ret;

  420. 

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

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

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

  424. 			nRows++;

  425. 		}

  426. 	}

  427. 	rewind(input);

  428. 

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

  430. 	double y[nRows*nColumns];

  431. 

  432. 	size_t i=0;

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

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

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

  436. 		}

  437. 		else{

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

  439. 			size_t j=0;

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

  441. 			j++;

  442. 			while(ret!=NULL){

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

  444. 				if(j<nColumns){

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

  446. 				}

  447. 				j++;

  448. 			}

  449. 			i++;

  450. 		}

  451. 	}

  452. 

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

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

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

  456. 		}

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

  458. 	}

  459. 	fclose(output);

  460. 

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

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

  463. 

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

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

  466. 	}

  467. 

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

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

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

  471. 	//}

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

  473. 	double dX0;

  474. 	int NN = nPts-2;

  475. 	

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

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

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

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

  480. 	}

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

  482. 

  483. 	gsl_interp_accel* acc;

  484. 	gsl_spline* spline;

  485. 	acc = gsl_interp_accel_alloc();

  486. 	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);

  487. 

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

  489. 

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

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

  492. 		}

  493. 

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

  495. 

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

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

  498. 		}

  499. 	}

  500. 	

  501. 	gsl_interp_accel_free(acc);

  502. 	gsl_spline_free(spline);

  503. 

  504. }

  505. 

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

  507. 	double mass=0.0e0;

  508. 	double rho;

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

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

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

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

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

  514. 	}

  515. 	return(mass);

  516. }

  517. 

  518. 

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

  520. 

  521. 	double rho,rhom;

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

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

  524. 	 * of kg. */

  525. 	psi(1)=ZERO;

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

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

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

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

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

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

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

  533. 		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;

  534. 	}

  535. 

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

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

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

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

  540. 	 * equations so that units match.*/

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

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

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

  544. 	}

  545. 	return(0);

  546. }

  547. 

  548. 

  549. int setInitialCondition(N_Vector* y, 

  550. 	    		N_Vector* ydot, 

  551.             		UserData data){

  552. 

  553. 	double* ydata;

  554. 	double* ydotdata;

  555. 	double* psidata;

  556. 	double* innerMassFractionsData, Rd, massDrop;

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

  558. 	double f=ZERO;

  559. 	double g=ZERO;

  560. 

  561. 	double perturb,rho;

  562. 	double epsilon=ZERO;

  563. 	int m,ier;

  564.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  565.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  566. 	innerMassFractionsData =  data->innerMassFractions;

  567. 	Rd = data->Rd;

  568. 	massDrop = data->massDrop;

  569. 

  570. 	//Mass of droplet

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

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

  573. 	if(data->adaptiveGrid){

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

  575. 	}

  576. 	else{

  577.   		psidata  = data->uniformGrid;

  578. 	}

  579. 

  580. 	m=data->metric;

  581. 

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

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

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

  585. 	}

  586. 

  587. 	//Define Grid:

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

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

  590. 	

  591.     //Initialize the R(i),T(i)(data->initialTemperature),Y(i,k),P(i)

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

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

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

  595. 		}else{

  596. 			if(i==1){

  597. 				R(i)=ZERO;

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

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

  600. 			}else{

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

  602. 			}

  603. 		}

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

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

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

  607. 		}

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

  609. 	}	

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

  611.     

  612. //    /********** test R(i) and the volumn between grids *****************/

  613. //    printf("Print the first 4 R(i)s and volumn between them: \n") ; 

  614. //    printf("Grids: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e \n",R(1),R(2),R(3),R(4));

  615. //    double v1,v2,v3,v4,v5; 

  616. //    v1 =pow(R(2),1+data->metric) - pow(R(1),1+data->metric);

  617. //    v2 =pow(R(3),1+data->metric) - pow(R(2),1+data->metric);

  618. //    v3 =pow(R(4),1+data->metric) - pow(R(3),1+data->metric);

  619. //    v4 =pow(R(5),1+data->metric) - pow(R(4),1+data->metric); 

  620. //    v5 =pow(R(6),1+data->metric) - pow(R(5),1+data->metric);

  621. //    printf("Volumn: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e,5th:%2.6e\n",v1,v2,v3,v4,v5) ; 

  622. 

  623. 

  624. 	double Tmax;

  625. 	double Tmin=data->initialTemperature;

  626. 	double w=data->mixingWidth;

  627. 	double YN2=ZERO;

  628. 	double YO2=ZERO;

  629. 	double YFuel,YOxidizer,sum;

  630. 

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

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

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

  634. 	//}

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

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

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

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

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

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

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

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

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

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

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

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

  647. 	//		}

  648. 	//	}

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

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

  651. 	//	}

  652. 	//}

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

  654. 	//	FILE* input;

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

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

  657. 	//		fclose(input);

  658. 	//	}

  659. 	//	else{

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

  661. 	//		return(-1);

  662. 	//	}

  663. 	//}

  664. 

  665. 	initializePsiGrid(ydata,psidata,data);

  666. 

  667. 	if(data->adaptiveGrid){

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

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

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

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

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

  673. 	//	}

  674. 	//	else{

  675. 	//	}

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

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

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

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

  680. 

  681. //            /**************** TEST THE data->grid->xOld *******************/

  682. //            double* ptr3 = data->grid->xOld ; 

  683. //            printf("SetInitialCondition function is called,before reGrid,Start print the first 5 elements of the xOld array : \n");

  684. //            printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr3[0],ptr3[1],ptr3[2],ptr3[3],ptr3[4]);

  685. 			

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

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

  688. 

  689. //            /**************** TEST THE data->grid->xOld *******************/

  690. //            double* ptr = data->grid->xOld ; 

  691. //            printf("SetInitialCondition function is called,after reGrid,Start print the first 5 elements of the xOld array : \n");

  692. //            printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr[0],ptr[1],ptr[2],ptr[3],ptr[4]);

  693. 

  694. 

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

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

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

  698. 		}

  699. 	}

  700. 	else{

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

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

  703. 		double psiNew[data->npts];

  704. 		

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

  706. 		psiNew[0] = 0;

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

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

  709. 		}

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

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

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

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

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

  715. 

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

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

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

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

  720. 		//}

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

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

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

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

  725. 	}

  726. 

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

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

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

  730. 	}

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

  732. 		/*Premixed Combustion: Equilibrium products comprise ignition

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

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

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

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

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

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

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

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

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

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

  743. 			}

  744. 		}

  745. 		if(data->dirichletOuter){

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

  747. 		}

  748. 	}

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

  750. 		FILE* input;

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

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

  753. 			fclose(input);

  754. 		}

  755. 		else{

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

  757. 			return(-1);

  758. 		}

  759. 		initializePsiGrid(ydata,psidata,data);

  760. 	}

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

  762. 		FILE* input;

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

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

  765. 			fclose(input);

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

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

  768. 			return(0);

  769. 		}

  770. 		else{

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

  772. 			return(-1);

  773. 		}

  774. 	}

  775. 

  776. 	if(data->reflectProblem){

  777. 		double temp;

  778. 		int j=1;

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

  780. 			temp=T(j);

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

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

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

  784. 				temp=Y(j,k);

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

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

  787. 			}

  788. 			j=j+1;

  789. 		}

  790. 	}

  791. 

  792. 	/*Floor small values to zero*/

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

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

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

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

  797. 			}

  798. 		}

  799. 	}

  800. 

  801. 	//Set grid to location of maximum curvature calculated by r instead of x:

  802. 	if(data->adaptiveGrid){

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

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

  805.         int maxCurvII = 0; 

  806.         maxCurvII = maxCurvIndexR(ydata,data->nt,data->nvar,data->nr,data->npts);

  807.         

  808.         data->grid->position = data->grid->x[maxCurvII];

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

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

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

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

  813.         

  814.         /******** Test the maxCurvPosition and related variables *******/

  815.         printf("The maxCurvPositition(based on r) is : %.6f,Temperature is : %.3f \n",data->grid->position,T(maxCurvII+1));

  816. 	}

  817. 

  818. 

  819. 	/*Ensure consistent boundary conditions*/	

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

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

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

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

  824. 	//}

  825. 

  826. 	return(0);

  827. }

  828. 

  829. inline double Qdot(double* t, 

  830. 	    	   double* x, 

  831. 	    	   double* ignTime, 

  832. 	    	   double* kernelSize, 

  833. 	    	   double* maxQdot){

  834. 	double qdot;

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

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

  837. 			qdot=(*maxQdot);

  838. 		}

  839. 		else{

  840. 			qdot=0.0e0;

  841. 		}

  842. 	}else{

  843. 		qdot=0.0e0;

  844. 	}

  845. 	return(qdot);

  846. }

  847. 

  848. inline void setGas(UserData data, 

  849. 	    double *ydata, 

  850. 	    size_t gridPoint){

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

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

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

  854. }

  855. 

  856. void getTransport(UserData data, 

  857. 		  double *ydata, 

  858. 		  size_t gridPoint,

  859. 		  double *rho,

  860. 		  double *lambda,

  861. 		  double YV[]){

  862. 

  863. 	double YAvg[data->nsp],

  864. 		 XLeft[data->nsp],

  865. 		 XRight[data->nsp],

  866. 		 gradX[data->nsp];

  867. 

  868. 	setGas(data,ydata,gridPoint);

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

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

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

  872.                                                                              

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

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

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

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

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

  878. 	}                                                                    

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

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

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

  882. 

  883. 

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

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

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

  887. 

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

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

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

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

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

  893. }

  894. 

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

  896. 

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

  898. 

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

  900.   	

  901. 	UserData data;

  902. 	data = (UserData)user_data;

  903. 	size_t npts=data->npts;

  904. 	size_t nsp=data->nsp;

  905. 	size_t k_bath = data->k_bath;

  906. 	size_t k_drop = data->k_drop;

  907. 

  908. 	ydata 	= N_VGetArrayPointer_OpenMP(y); 

  909. 	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 

  910. 	resdata = N_VGetArrayPointer_OpenMP(res);

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

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

  913. 	}else{

  914. 		psidata = data->uniformGrid;

  915. 	}

  916. 

  917. 	innerMassFractionsData =  data->innerMassFractions;

  918. 

  919. 	/* Grid stencil:*/

  920. 

  921. 	/*-------|---------*---------|---------*---------|-------*/

  922. 	/*-------|---------*---------|---------*---------|-------*/

  923. 	/*-------|---------*---------|---------*---------|-------*/

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

  925. 	/*-------|---------*---------|---------*---------|-------*/

  926. 	/*-------|---------*---------|---------*---------|-------*/

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

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

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

  930. 

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

  932. 	 * calculated values:

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

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

  935. 	 * m 		: exponent that determines geometry;

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

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

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

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

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

  941. 	 * Tav		: average temperature between two points.

  942. 	 * Pav		: average pressure between two points.

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

  944. 	 * Xgradhalf	: mole fraction gradient at j.

  945. 	 * Cpb		: mass based bulk specific heat.

  946. 	 * tranTerm	: transient terms.

  947. 	 * advTerm	: advection terms.

  948. 	 * diffTerm	: diffusion terms.

  949. 	 * srcTerm	: source terms.

  950. 	 */

  951. 

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

  953. 		 rho, 

  954. 		 rhophalf, lambdaphalf, YVphalf[nsp],

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

  956. 		 tranTerm, diffTerm,    srcTerm, advTerm,

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

  958. 	

  959. 	/*Aliases for difference coefficients:*/

  960. 	double cendfm, cendfc, cendfp;

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

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

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

  964. 	//double mass, mdotIn;

  965. 	double mass, massDrop;

  966. 	double sum, sum1, sum2, sum3;

  967. 

  968. 	size_t j,k;

  969. 	int m;

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

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

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

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

  974. 

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

  976. 	setGas(data,ydata,1);

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

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

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

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

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

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

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

  984.     

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

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

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

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

  989. 

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

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

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

  993. 	

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

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

  996.     

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

  998.     //rhol = 997.0;	

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

  1000.     

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

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

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

  1004. 

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

  1006.     

  1007. 

  1008. 	/*******************************************************************/

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

  1010. 

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

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

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

  1014. 	areamhalfsq= areamhalf*areamhalf;                                   

  1015.                                                                              

  1016. 	/*******************************************************************/

  1017. 

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

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

  1020. 

  1021. 

  1022. 	/*Mass:*/                                                            

  1023. 	if(data->quasiSteady){

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

  1025. 	}else{

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

  1027. 	}

  1028.                                                                         

  1029. 	/*Species:*/                                                    

  1030. 	sum=ZERO;

  1031. 	

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

  1033. 	//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) 

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

  1035. 	//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) 

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

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

  1038. 	

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

  1040.     //Antoine parameter from NIST website

  1041.     double p_i=0.0 ; 

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

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

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

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

  1046. 	

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

  1048. 

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

  1050. 

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

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

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

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

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

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

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

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

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

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

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

  1062. 			//}

  1063. 			//else{

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

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

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

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

  1068. 			//	 * fix this?*/

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

  1070. 			//}

  1071. 		}

  1072. 	}                                                                 

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

  1074. 	

  1075. 	/*Energy:*/                                                     

  1076. 	if(data->dirichletInner){

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

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

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

  1080. 	}else{

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

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

  1083. 		//	/ (massDrop * Cpl);

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

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

  1086. 			/ (massDrop * Cpl);

  1087. 

  1088. 

  1089. 

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

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

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

  1093. 	}

  1094. 

  1095. 	/*Pressure:*/                                                        

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

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

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

  1099.                                                                              

  1100. 		/*evaluate various mesh differences*///                      

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

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

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

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

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

  1106.                                                                              

  1107. 		/*evaluate various central difference coefficients*/       

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

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

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

  1111. 		/**************************************************/

  1112.                                                                    

  1113.                                                                   

  1114. 		/*evaluate properties at j*************************/

  1115. 		setGas(data,ydata,j);

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

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

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

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

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

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

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

  1123.                                                                            

  1124. 		/*evaluate properties at p*************************/

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

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

  1127. 		areaphalfsq= areaphalf*areaphalf;		    

  1128. 		/**************************************************///     

  1129. 

  1130. 		/*Evaporating Mass*/

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

  1132. 

  1133. 		/*Mass:*/

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

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

  1136. 

  1137. 		/*Energy:*/

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

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

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

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

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

  1143. 		/*Notes: 

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

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

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

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

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

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

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

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

  1152. 		 * units.

  1153. 		 * */

  1154. 

  1155. 		//enthalpy formulation:

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

  1157. 		tranTerm = Tdot(j);

  1158. 		sum=ZERO;

  1159. 		sum1=ZERO;

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

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

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

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

  1164. 		}

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

  1166. 		sum1=sum1*Cantera::GasConstant;

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

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

  1169. 			     /(dpsiav*Cpb) )

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

  1171. 			              +cendfc*T(j)

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

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

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

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

  1176. 			+advTerm 

  1177. 			-diffTerm

  1178. 			+srcTerm;

  1179. 

  1180. 	//	//energy formulation:

  1181. 	//	tranTerm = Tdot(j);

  1182. 	//	sum=ZERO;

  1183. 	//	sum1=ZERO;

  1184. 	//	sum2=ZERO;

  1185. 	//	sum3=ZERO;

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

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

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

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

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

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

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

  1193. 	//	}

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

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

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

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

  1198. 	//		     /(dpsiav*Cvb) )

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

  1200. 	//		              +cendfc*T(j)

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

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

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

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

  1205. 	//	Tres(j)= tranTerm

  1206. 	//		+advTerm 

  1207. 	//		-diffTerm

  1208. 	//		+srcTerm;

  1209. 

  1210. 		/*Species:*/

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

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

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

  1214. 

  1215. 		sum=ZERO;

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

  1217. 			if(k!=k_bath){

  1218. 				tranTerm  = Ydot(j,k);

  1219. 				diffTerm  = (YVphalf(k)*areaphalf

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

  1221. 				srcTerm   = wdot(k)

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

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

  1224. 

  1225. 				Yres(j,k)= tranTerm 

  1226. 					  +advTerm

  1227. 					  +diffTerm

  1228. 					  -srcTerm;

  1229. 

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

  1231. 			}

  1232. 		}

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

  1234. 

  1235. 		/*Pressure:*/

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

  1237. 

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

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

  1240. 		areamhalf=areaphalf;                                

  1241. 		areamhalfsq=areaphalfsq;                            

  1242. 		aream=area;

  1243. 		rhom=rho;

  1244. 		rhomhalf=rhophalf;                                  

  1245. 		lambdamhalf=lambdaphalf;                            

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

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

  1248. 		}                                                   

  1249. 		/**************************************************/

  1250. 		                                                    

  1251. 	}                                                           

  1252. 	/*******************************************************************///

  1253. 

  1254. 

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

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

  1257. 	setGas(data,ydata,npts);

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

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

  1260. 

  1261. 	/*Mass:*/                                                          

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

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

  1264. 

  1265. 	/*Energy:*/                                             

  1266. 	if(data->dirichletOuter){

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

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

  1269. 	}

  1270. 	else{

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

  1272. 	}

  1273.                                                                 

  1274. 	/*Species:*/                                            

  1275. 	sum=ZERO;

  1276. 	if(data->dirichletOuter){

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

  1278. 			if(k!=k_bath){

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

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

  1281. 			}

  1282. 		}                                                       

  1283. 	}

  1284. 	else{

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

  1286. 			if(k!=k_bath){

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

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

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

  1290. 			}

  1291. 		}                                                       

  1292. 	}

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

  1294.                                                                 

  1295. 	

  1296. 	/*Pressure:*/                                                      

  1297. 	if(data->constantPressure){

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

  1299. 	}

  1300. 	else{

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

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

  1303. 	}

  1304. 

  1305. 	/*Evaporating Mass*/

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

  1307. 

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

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

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

  1311. 	//	//}

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

  1313. 	//}

  1314. 

  1315. 	return(0);

  1316. 

  1317. }

  1318. 

  1319. void printSpaceTimeHeader(UserData data)

  1320. {

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

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

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

  1324. 	}

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

  1326. 

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

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

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

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

  1331. 	}

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

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

  1334. }

  1335. 

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

  1337. {

  1338. 	double *ydata,*psidata;

  1339.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1340. 

  1341. 	if(data->adaptiveGrid){

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

  1343. 	}else{

  1344.   		psidata  = data->uniformGrid;

  1345. 	}

  1346. 

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

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

  1349. 		if(i==0){

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

  1351. 		}

  1352. 		else{

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

  1354. 		}

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

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

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

  1358. 			}else{

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

  1360. 			}

  1361. 		}

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

  1363. 	}

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

  1365. }

  1366. 

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

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

  1369. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1370. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1371. 	if(data->adaptiveGrid){

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

  1373. 	}else{

  1374.   		psidata  = data->uniformGrid;

  1375. 	}

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

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

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

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

  1380. }

  1381. 

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

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

  1384. 	double t;

  1385. 	if(data->adaptiveGrid){

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

  1387. 	}else{

  1388.   		psidata  = data->uniformGrid;

  1389. 	}

  1390. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1391. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

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

  1393. 	data->tNow=t;

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

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

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

  1397. 	if(data->adaptiveGrid){

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

  1399. 	}

  1400. }

  1401. 

  1402. void printGlobalHeader(UserData data)

  1403. {

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

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

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

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

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

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

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

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

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

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

  1414. }

  1415. //

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

  1417. //{

  1418. //	double *ydotdata,*psidata;

  1419. //  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);

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

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

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

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

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

  1425. //		}

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

  1427. //	}

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

  1429. //}

  1430. //

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

  1432. {

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

  1434. 	innerMassFractionsData =  data->innerMassFractions;

  1435. 

  1436. 	if(data->adaptiveGrid){

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

  1438. 	}else{

  1439.   		psidata  = data->uniformGrid;

  1440. 	}

  1441. 

  1442. 	double TAvg, RAvg, YAvg, psiAvg;

  1443.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1444.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1445. 	TAvg=data->isotherm;

  1446. 	double sum=ZERO;

  1447. 	double dpsim,area,aream,drdt;

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

  1449. 

  1450. 

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

  1452. 	size_t j=1;

  1453. 	size_t jj=1;

  1454. 	size_t jjj=1;

  1455. 	double wdot[data->nsp];

  1456. 	double wdotMax=0.0e0;

  1457. 	double advTerm=0.0e0;

  1458. 	psiAvg=0.0e0;

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

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

  1461. 			j=j+1;

  1462. 		}

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

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

  1465. 			jj=jj+1;

  1466. 		}

  1467. 

  1468. 		setGas(data,ydata,jj);

  1469. 		Tu=T(jj);

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

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

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

  1473. 		gamma=Cpb/Cvb;

  1474. 		

  1475. 	}

  1476. 	else{

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

  1478. 			j=j+1;

  1479. 		}

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

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

  1482. 			jj=jj+1;

  1483. 		}

  1484. 

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

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

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

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

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

  1490. 		gamma=Cpb/Cvb;

  1491. 	}

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

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

  1494. 	}

  1495. 	else{

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

  1497. 	}

  1498. 

  1499. 	////Experimental burning speed calculation:	

  1500. 	//int nMax=0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1519. 	////		}

  1520. 	////	}

  1521. 	////}

  1522. 	////else{

  1523. 	////	psiAvg=0.0e0;

  1524. 	////}

  1525. 

  1526. 

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

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

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

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

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

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

  1533. 

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

  1535. 	 * the consumption of O2*/

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

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

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

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

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

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

  1542. 		aream=area;

  1543. 	}

  1544. 

  1545. 	//double maxOH,maxHO2;

  1546. 	//maxOH=0.0e0;

  1547. 	//maxHO2=0.0e0;

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

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

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

  1551. 	//	}

  1552. 	//}

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

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

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

  1556. 	//	}

  1557. 	//}

  1558. 

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

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

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

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

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

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

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

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

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

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

  1569. }

  1570. 

  1571. //

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

  1573. //{

  1574. //	double *ydata,*psidata;

  1575. //  	ydata    = N_VGetArrayPointer_OpenMP(y);

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

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

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

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

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

  1581. //		}

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

  1583. //	}

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

  1585. //}

  1586. //

  1587. //

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

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

  1590. ////  	double *ydata;

  1591. ////  	double *ydotdata;

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

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

  1594. ////

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

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

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

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

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

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

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

  1602. ////

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

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

  1605. ////	}

  1606. ////}

  1607. 

  1608. 

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

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

  1611. 	size_t i, k, nsp,npts ; 

  1612. 	nsp = data->nsp ; 

  1613. 	npts = data->npts ; 

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

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

  1616. 

  1617. 	double* ydata ; 

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

  1619. 	

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

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

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

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

  1624. 		

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

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

  1627. 		}

  1628. 

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

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

  1631. 		}

  1632. 

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

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

  1635. 		}

  1636. 

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

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

  1639. 		}

  1640.         

  1641. 

  1642. 

  1643. 	}

  1644.     

  1645. }

  1646. 

  1647. 

  1648. void printTimescaleHeader(UserData data)

  1649. {

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

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

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

  1653. 	}

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

  1655. 

  1656. 

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

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

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

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

  1661. 	}

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

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

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

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

  1666. }

  1667. 

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

  1669. {

  1670. 	double* ydata; 

  1671. 	ydata = N_VGetArrayPointer_OpenMP(*y); 

  1672. 

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

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

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

  1676.         

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

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

  1679. 		}

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

  1681. 	}

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

  1683. 

  1684. }

  1685. 

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

  1687. {

  1688.     double* ydata;

  1689.     ydata = N_VGetArrayPointer_OpenMP(*y);

  1690.     //double sum = 0.00; 

  1691.     size_t k_bath = data->k_bath;

  1692.    

  1693. 	/*Floor small values to zero*/

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

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

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

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

  1698. 			}

  1699. 		}

  1700. 	}

  1701.     

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

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

  1704.         double sum = 0.00;

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

  1706. 			if(k!=k_bath){

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

  1708. 			}

  1709. 		}

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

  1711. 	}

  1712.     

  1713. 

  1714. }

  1715. 

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

  1717. {

  1718.     double* ydata;

  1719.     realtype* atolvdata; 

  1720.     ydata = N_VGetArrayPointer_OpenMP(*y);

  1721.     atolvdata = N_VGetArrayPointer_OpenMP(*atolv);

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

  1723.     double maxT=0.00, ignT=1800;

  1724.     relTol = 1.0e-03 ; 

  1725.     radTol = 1.0e-05 ; 

  1726.     tempTol = 1.0e-01 ; 

  1727.     presTol = 1.0e00 ; 

  1728.     massFracTol = 1.0e-08; 

  1729.     bathGasTol = 1.0e-04 ; 

  1730.     mdotTol = 1.0e-09 ; 

  1731. 

  1732. 

  1733.     /*Get the maximum Tempture*/

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

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

  1736.             maxT = T(i);

  1737.         }

  1738.     }

  1739. 

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

  1741.     if(maxT >= ignT){

  1742.         data->relativeTolerance = relTol;

  1743. 

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

  1745.             atolT(i) = tempTol; 

  1746.             atolR(i) = radTol ; 

  1747.             atolP(i) = presTol ;

  1748.             atolMdot(i) = mdotTol ; 

  1749.             

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

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

  1752.                     atolY(i,k) = massFracTol;

  1753.                 }else{

  1754.                     atolY(i,k) = bathGasTol;

  1755.                 }

  1756.             }

  1757.         }

  1758.     }

  1759. 

  1760. }