USC
/
LTORC
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
A low Mach, 1D, reacting flow code.
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.
1 Commit
1 Branch
58KB
Tree: 48cca4770f
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from '48cca4770f'
${ noResults }
LTORC/residue.cpp

1304 lines
37KB
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. void updateSolution(double* y, double* ydot, const size_t nvar,

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

  114. 

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

  116. 

  117. 	gsl_interp_accel* acc;

  118. 	gsl_spline* spline;

  119. 	acc = gsl_interp_accel_alloc();

  120. 	spline = gsl_spline_alloc(gsl_interp_steffen, nPts);

  121. 

  122. 	gsl_interp_accel* accdot;

  123. 	gsl_spline* splinedot;

  124. 	accdot    = gsl_interp_accel_alloc();

  125. 	splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);

  126. 

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

  128. 

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

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

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

  132. 		}

  133. 

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

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

  136. 

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

  138. 

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

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

  141. 		}

  142. 	}

  143. 

  144. 	//Exploring "fixing" boundary conditions:

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

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

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

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

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

  150. 	//}

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

  152. 	

  153. 	gsl_interp_accel_free(acc);

  154. 	gsl_spline_free(spline);

  155. 

  156. 	gsl_interp_accel_free(accdot);

  157. 	gsl_spline_free(splinedot);

  158. }

  159. 

  160. //Locate bath gas:

  161. size_t BathGasIndex(UserData data){

  162. 	size_t index=0;

  163. 	double max1;

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

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

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

  167. 		if(max1>=max){

  168. 			max=max1;

  169. 			index=k;

  170. 		}

  171. 	}

  172. 	return(index+1);

  173. }

  174. 

  175. //Locate Oxidizer:

  176. size_t oxidizerIndex(UserData data){

  177. 	size_t index=0;

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

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

  180. 			index=k;

  181. 		}

  182. 	}

  183. 	return(index);

  184. }

  185. 

  186. //Locate OH:

  187. size_t OHIndex(UserData data){

  188. 	size_t index=0;

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

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

  191. 			index=k;

  192. 		}

  193. 	}

  194. 	return(index);

  195. }

  196. 

  197. //Locate HO2:

  198. size_t HO2Index(UserData data){

  199. 	size_t index=0;

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

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

  202. 			index=k;

  203. 		}

  204. 	}

  205. 	return(index);

  206. }

  207. 

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

  209. 	double *iddata;

  210.   	N_VConst(ONE, *id);

  211.   	iddata   = N_VGetArrayPointer_OpenMP(*id);

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

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

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

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

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

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

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

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

  220. 		Rid(i)=ZERO;

  221. 		Pid(i)=ZERO;

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

  223. 	}	

  224. 	Rid(1)=ONE;

  225. 	Tid(1)=ZERO;

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

  227. 	if(data->constantPressure){

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

  229. 	}

  230. 	else{

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

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

  233. 	}

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

  235. 		Yid(1,k)=ZERO;

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

  237. 	}

  238. 	return(0);

  239. }

  240. 

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

  242. 	switch (*i) {

  243. 		case 0:

  244. 			return(ONE);

  245. 		case 1:

  246. 			return(x);

  247. 		case 2:

  248. 			return(x*x);

  249. 		default:

  250. 			return(ONE);

  251. 	}

  252. }

  253. 

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

  255. 

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

  257. 

  258. 	size_t bufLen=10000;

  259. 	size_t nRows=0;

  260. 	size_t nColumns=nvar;

  261. 

  262. 	char buf[bufLen];

  263. 	char buf1[bufLen];

  264. 	char comment[1];

  265. 	char *ret;

  266. 

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

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

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

  270. 			nRows++;

  271. 		}

  272. 	}

  273. 	rewind(input);

  274. 

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

  276. 	double y[nRows*nColumns];

  277. 

  278. 	size_t i=0;

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

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

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

  282. 		}

  283. 		else{

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

  285. 			size_t j=0;

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

  287. 			j++;

  288. 			while(ret!=NULL){

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

  290. 				if(j<nColumns){

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

  292. 				}

  293. 				j++;

  294. 			}

  295. 			i++;

  296. 		}

  297. 	}

  298. 

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

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

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

  302. 		}

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

  304. 	}

  305. 	fclose(output);

  306. 

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

  308. 

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

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

  311. 	}

  312. 

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

  314. 	for (size_t j = 0; j < nPts; j++) {

  315. 		xNew[j]=(double)j*dx;

  316. 	}

  317. 

  318. 	gsl_interp_accel* acc;

  319. 	gsl_spline* spline;

  320. 	acc = gsl_interp_accel_alloc();

  321. 	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);

  322. 

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

  324. 

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

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

  327. 		}

  328. 

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

  330. 

  331. 		for (size_t k = 1; k < nPts-1; k++) {

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

  333. 		}

  334. 	}

  335. 	

  336. 	gsl_interp_accel_free(acc);

  337. 	gsl_spline_free(spline);

  338. 

  339. }

  340. 

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

  342. 	double mass=0.0e0;

  343. 	double rho;

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

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

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

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

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

  349. 	}

  350. 	return(mass);

  351. }

  352. 

  353. 

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

  355. 

  356. 	double rho;

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

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

  359. 	 * of kg. */

  360. 	psi(1)=ZERO;

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

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

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

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

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

  366. 	}

  367. 

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

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

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

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

  372. 	 * equations so that units match.*/

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

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

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

  376. 	}

  377. 	return(0);

  378. }

  379. 

  380. 

  381. int setInitialCondition(N_Vector* y, 

  382. 	    		N_Vector* ydot, 

  383.             		UserData data){

  384. 

  385. 	double* ydata;

  386. 	double* ydotdata;

  387. 	double* psidata;

  388. 	double* innerMassFractionsData;

  389. 	double f=ZERO;

  390. 	double g=ZERO;

  391. 

  392. 	double perturb,rho;

  393. 	double epsilon=ZERO;

  394. 	int m,ier;

  395.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  396.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  397. 	innerMassFractionsData =  data->innerMassFractions;

  398. 

  399. 	if(data->adaptiveGrid){

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

  401. 	}

  402. 	else{

  403.   		psidata  = data->uniformGrid;

  404. 	}

  405. 

  406. 	m=data->metric;

  407. 

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

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

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

  411. 	}

  412. 

  413. 	//Define Grid:

  414. 	//double Rmin=1e-03*data->domainLength;

  415. 	double Rmin=0.0e0;

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

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

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

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

  420. 			R(i)=Rmin+(double)((i-1)*dR);

  421. 		}else{

  422. 			if(i==1){

  423. 				R(i)=ZERO;

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

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

  426. 			}else{

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

  428. 			}

  429. 		}

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

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

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

  433. 		}

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

  435. 	}	

  436. 	R(data->npts)=data->domainLength;

  437. 

  438. 	double Tmax;

  439. 	double Tmin=data->initialTemperature;

  440. 	double w=data->mixingWidth;

  441. 	double YN2=ZERO;

  442. 	double YO2=ZERO;

  443. 	double YFuel,YOxidizer,sum;

  444. 

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

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

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

  448. 	//}

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

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

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

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

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

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

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

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

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

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

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

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

  461. 	//		}

  462. 	//	}

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

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

  465. 	//	}

  466. 	//}

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

  468. 	//	FILE* input;

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

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

  471. 	//		fclose(input);

  472. 	//	}

  473. 	//	else{

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

  475. 	//		return(-1);

  476. 	//	}

  477. 	//}

  478. 

  479. 	initializePsiGrid(ydata,psidata,data);

  480. 

  481. 	if(data->adaptiveGrid){

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

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

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

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

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

  487. 	//	}

  488. 	//	else{

  489. 	//	}

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

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

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

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

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

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

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

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

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

  499. 		}

  500. 	}

  501. 	else{

  502. 		double psiNew[data->npts];

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

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

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

  506. 		}

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

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

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

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

  511. 	}

  512. 

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

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

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

  516. 	}

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

  518. 		/*Premixed Combustion: Equilibrium products comprise ignition

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

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

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

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

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

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

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

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

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

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

  529. 			}

  530. 		}

  531. 		if(data->dirichletOuter){

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

  533. 		}

  534. 	}

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

  536. 		FILE* input;

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

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

  539. 			fclose(input);

  540. 		}

  541. 		else{

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

  543. 			return(-1);

  544. 		}

  545. 	}

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

  547. 		FILE* input;

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

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

  550. 			fclose(input);

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

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

  553. 			return(0);

  554. 		}

  555. 		else{

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

  557. 			return(-1);

  558. 		}

  559. 	}

  560. 

  561. 	if(data->reflectProblem){

  562. 		double temp;

  563. 		int j=1;

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

  565. 			temp=T(j);

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

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

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

  569. 				temp=Y(j,k);

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

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

  572. 			}

  573. 			j=j+1;

  574. 		}

  575. 	}

  576. 

  577. 	/*Floor small values to zero*/

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

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

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

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

  582. 			}

  583. 		}

  584. 	}

  585. 

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

  587. 	if(data->adaptiveGrid){

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

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

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

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

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

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

  594. 	}

  595. 

  596. 	/*Ensure consistent boundary conditions*/	

  597. 	T(1)=T(2);

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

  599. 		Y(1,k)=Y(2,k);

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

  601. 	}

  602. 

  603. 

  604. 	return(0);

  605. }

  606. 

  607. inline double Qdot(double* t, 

  608. 	    	   double* x, 

  609. 	    	   double* ignTime, 

  610. 	    	   double* kernelSize, 

  611. 	    	   double* maxQdot){

  612. 	double qdot;

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

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

  615. 			qdot=(*maxQdot);

  616. 		}

  617. 		else{

  618. 			qdot=0.0e0;

  619. 		}

  620. 	}else{

  621. 		qdot=0.0e0;

  622. 	}

  623. 	return(qdot);

  624. }

  625. 

  626. inline void setGas(UserData data, 

  627. 	    double *ydata, 

  628. 	    size_t gridPoint){

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

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

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

  632. }

  633. 

  634. void getTransport(UserData data, 

  635. 		  double *ydata, 

  636. 		  size_t gridPoint,

  637. 		  double *rho,

  638. 		  double *lambda,

  639. 		  double YV[]){

  640. 

  641. 	double YAvg[data->nsp],

  642. 		 XLeft[data->nsp],

  643. 		 XRight[data->nsp],

  644. 		 gradX[data->nsp];

  645. 

  646. 	setGas(data,ydata,gridPoint);

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

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

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

  650.                                                                              

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

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

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

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

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

  656. 	}                                                                    

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

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

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

  660. 

  661. 

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

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

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

  665. 

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

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

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

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

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

  671. }

  672. 

  673. int residue(double t, 

  674. 	    N_Vector y, 

  675. 	    N_Vector ydot, 

  676. 	    N_Vector res, 

  677.             void *user_data){

  678. 

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

  680. 

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

  682.   	

  683. 	UserData data;

  684. 	data = (UserData)user_data;

  685. 	size_t npts=data->npts;

  686. 	size_t nsp=data->nsp;

  687. 	size_t k_bath = data->k_bath;

  688. 

  689. 	ydata 	= N_VGetArrayPointer_OpenMP(y); 

  690. 	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 

  691. 	resdata = N_VGetArrayPointer_OpenMP(res);

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

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

  694. 	}else{

  695. 		psidata = data->uniformGrid;

  696. 	}

  697. 

  698. 	innerMassFractionsData =  data->innerMassFractions;

  699. 

  700. 	/* Grid stencil:*/

  701. 

  702. 	/*-------|---------*---------|---------*---------|-------*/

  703. 	/*-------|---------*---------|---------*---------|-------*/

  704. 	/*-------|---------*---------|---------*---------|-------*/

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

  706. 	/*-------|---------*---------|---------*---------|-------*/

  707. 	/*-------|---------*---------|---------*---------|-------*/

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

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

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

  711. 

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

  713. 	 * calculated values:

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

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

  716. 	 * m 		: exponent that determines geometry;

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

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

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

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

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

  722. 	 * Tav		: average temperature between two points.

  723. 	 * Pav		: average pressure between two points.

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

  725. 	 * Xgradhalf	: mole fraction gradient at j.

  726. 	 * Cpb		: mass based bulk specific heat.

  727. 	 * tranTerm	: transient terms.

  728. 	 * advTerm	: advection terms.

  729. 	 * diffTerm	: diffusion terms.

  730. 	 * srcTerm	: source terms.

  731. 	 */

  732. 

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

  734. 		 rho, 

  735. 		 rhophalf, lambdaphalf, YVphalf[nsp],

  736. 		 Cpb, Cvb,  Cp[nsp],   wdot[nsp],   enthalpy[nsp], energy[nsp],

  737. 		 tranTerm, diffTerm,    srcTerm, advTerm,

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

  739. 	

  740. 	/*Aliases for difference coefficients:*/

  741. 	double cendfm, cendfc, cendfp;

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

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

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

  745. 	double mass, mdotIn;

  746. 	double sum, sum1, sum2, sum3;

  747. 

  748. 	size_t j,k;

  749. 	int m;

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

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

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

  753. 

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

  755. 	setGas(data,ydata,1);

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

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

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

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

  760. 

  761. 	/*******************************************************************/

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

  763. 

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

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

  766. 	areamhalfsq= areamhalf*areamhalf;                                   

  767.                                                                              

  768. 	/*******************************************************************/

  769. 

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

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

  772. 

  773. 	/*Mass:*/                                                            

  774. 	Rres(1)=Rdot(1);

  775. 

  776. 	/*Energy:*/                                                     

  777. 

  778. 	if(data->dirichletInner){

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

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

  781. 	}

  782. 	else{

  783. 		Tres(1)=T(2)-T(1);

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

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

  786. 	}

  787.                                                                         

  788. 	/*Species:*/                                                    

  789. 	sum=ZERO;

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

  791. 		if(k!=k_bath){

  792. 			if(fabs(mdotIn)>1e-14){

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

  794. 					  Y(1,k)-

  795. 					  (YVmhalf(k)*areamhalf)/mdotIn;

  796. 			}

  797. 			else{

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

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

  800. 				/*The below flux boundary condition makes the

  801. 				 * problem more prone to diverge. How does one

  802. 				 * fix this?*/

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

  804. 			}

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

  806. 		}

  807. 	}                                                                 

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

  809. 

  810.                                                                              

  811. 	/*Pressure:*/                                                        

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

  813. 

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

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

  816.                                                                              

  817. 		/*evaluate various mesh differences*///                      

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

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

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

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

  822. 		/***********************************///                      

  823.                                                                              

  824. 		/*evaluate various central difference coefficients*/       

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

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

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

  828. 		/**************************************************/

  829.                                                                    

  830.                                                                   

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

  832. 		setGas(data,ydata,j);

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

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

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

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

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

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

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

  840.                                                                            

  841. 		/*evaluate properties at p*************************/

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

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

  844. 		areaphalfsq= areaphalf*areaphalf;		    

  845. 		/**************************************************///     

  846. 

  847. 		/*Mass:*/

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

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

  850. 

  851. 		/*Energy:*/

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

  853. 		 * 	   + (1/cₚ)(∂/∂ψ)(λρA²∂T/∂ψ) 

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

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

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

  857. 		/*Notes: 

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

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

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

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

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

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

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

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

  866. 		 * units.

  867. 		 * */

  868. 

  869. 		//enthalpy formulation:

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

  871. 		sum=ZERO;

  872. 		sum1=ZERO;

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

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

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

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

  877. 		}

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

  879. 		sum1=sum1*Cantera::GasConstant;

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

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

  882. 			     /(dpsiav*Cpb) )

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

  884. 			              +cendfc*T(j)

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

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

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

  888. 		Tres(j)= tranTerm

  889. 			+advTerm 

  890. 			-diffTerm

  891. 			+srcTerm;

  892. 

  893. 	//	//energy formulation:

  894. 	//	tranTerm = Tdot(j);

  895. 	//	sum=ZERO;

  896. 	//	sum1=ZERO;

  897. 	//	sum2=ZERO;

  898. 	//	sum3=ZERO;

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

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

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

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

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

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

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

  906. 	//	}

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

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

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

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

  911. 	//		     /(dpsiav*Cvb) )

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

  913. 	//		              +cendfc*T(j)

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

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

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

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

  918. 	//	Tres(j)= tranTerm

  919. 	//		+advTerm 

  920. 	//		-diffTerm

  921. 	//		+srcTerm;

  922. 

  923. 		/*Species:*/

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

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

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

  927. 

  928. 		sum=ZERO;

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

  930. 			if(k!=k_bath){

  931. 				tranTerm  = Ydot(j,k);

  932. 				diffTerm  = (YVphalf(k)*areaphalf

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

  934. 				srcTerm   = wdot(k)

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

  936. 				advTerm   = (mdotIn*(Y(j,k)-Y(j-1,k))/dpsim);

  937. 

  938. 				Yres(j,k)= tranTerm 

  939. 					  +advTerm

  940. 					  +diffTerm

  941. 					  -srcTerm;

  942. 

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

  944. 			}

  945. 		}

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

  947. 

  948. 		/*Pressure:*/

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

  950. 

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

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

  953. 		areamhalf=areaphalf;                                

  954. 		areamhalfsq=areaphalfsq;                            

  955. 		aream=area;

  956. 		rhom=rho;

  957. 		rhomhalf=rhophalf;                                  

  958. 		lambdamhalf=lambdaphalf;                            

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

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

  961. 		}                                                   

  962. 		/**************************************************/

  963. 		                                                    

  964. 	}                                                           

  965. 	/*******************************************************************///

  966. 

  967. 

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

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

  970. 	setGas(data,ydata,npts);

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

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

  973. 

  974. 	/*Mass:*/                                                          

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

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

  977. 

  978. 	/*Energy:*/                                             

  979. 	if(data->dirichletOuter){

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

  981. 	}

  982. 	else{

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

  984. 	}

  985.                                                                 

  986. 	/*Species:*/                                            

  987. 	sum=ZERO;

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

  989. 		if(k!=k_bath){

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

  991. 		    	Yres(npts,k)=YVmhalf(k);         

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

  993. 		}

  994. 	}                                                       

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

  996.                                                                 

  997. 	

  998. 	/*Pressure:*/                                                      

  999. 	if(data->constantPressure){

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

  1001. 	}

  1002. 	else{

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

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

  1005. 	}

  1006. 

  1007. 

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

  1009. 

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

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

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

  1013. 	//	//}

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

  1015. 	//}

  1016. 

  1017. 	return(0);

  1018. 

  1019. }

  1020. 

  1021. void printSpaceTimeHeader(UserData data)

  1022. {

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

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

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

  1026. 	}

  1027. 	fprintf((data->output), "\n");

  1028. 

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

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

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

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

  1033. 	}

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

  1035. }

  1036. 

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

  1038. {

  1039. 	double *ydata,*psidata;

  1040.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1041. 

  1042. 	if(data->adaptiveGrid){

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

  1044. 	}else{

  1045.   		psidata  = data->uniformGrid;

  1046. 	}

  1047. 

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

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

  1050. 		if(i==0){

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

  1052. 		}

  1053. 		else{

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

  1055. 		}

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

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

  1058. 		}

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

  1060. 	}

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

  1062. }

  1063. 

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

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

  1066. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1067. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1068. 	if(data->adaptiveGrid){

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

  1070. 	}else{

  1071.   		psidata  = data->uniformGrid;

  1072. 	}

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

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

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

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

  1077. }

  1078. 

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

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

  1081. 	double t;

  1082. 	if(data->adaptiveGrid){

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

  1084. 	}else{

  1085.   		psidata  = data->uniformGrid;

  1086. 	}

  1087. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1088. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

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

  1090. 	data->tNow=t;

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

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

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

  1094. 	if(data->adaptiveGrid){

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

  1096. 	}

  1097. }

  1098. 

  1099. void printGlobalHeader(UserData data)

  1100. {

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

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

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

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

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

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

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

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

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

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

  1111. }

  1112. //

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

  1114. //{

  1115. //	double *ydotdata,*psidata;

  1116. //  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);

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

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

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

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

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

  1122. //		}

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

  1124. //	}

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

  1126. //}

  1127. //

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

  1129. {

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

  1131. 	innerMassFractionsData =  data->innerMassFractions;

  1132. 

  1133. 	if(data->adaptiveGrid){

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

  1135. 	}else{

  1136.   		psidata  = data->uniformGrid;

  1137. 	}

  1138. 

  1139. 	double TAvg, RAvg, YAvg, psiAvg;

  1140.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1141.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1142. 	TAvg=data->isotherm;

  1143. 	double sum=ZERO;

  1144. 	double dpsim,area,aream,drdt;

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

  1146. 

  1147. 

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

  1149. 	size_t j=1;

  1150. 	size_t jj=1;

  1151. 	size_t jjj=1;

  1152. 	double wdot[data->nsp];

  1153. 	double wdotMax=0.0e0;

  1154. 	double advTerm=0.0e0;

  1155. 	psiAvg=0.0e0;

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

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

  1158. 			j=j+1;

  1159. 		}

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

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

  1162. 			jj=jj+1;

  1163. 		}

  1164. 

  1165. 		setGas(data,ydata,jj);

  1166. 		Tu=T(jj);

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

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

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

  1170. 		gamma=Cpb/Cvb;

  1171. 		

  1172. 	}

  1173. 	else{

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

  1175. 			j=j+1;

  1176. 		}

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

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

  1179. 			jj=jj+1;

  1180. 		}

  1181. 

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

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

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

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

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

  1187. 		gamma=Cpb/Cvb;

  1188. 	}

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

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

  1191. 	}

  1192. 	else{

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

  1194. 	}

  1195. 

  1196. 	////Experimental burning speed calculation:	

  1197. 	//int nMax=0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1216. 	////		}

  1217. 	////	}

  1218. 	////}

  1219. 	////else{

  1220. 	////	psiAvg=0.0e0;

  1221. 	////}

  1222. 

  1223. 

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

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

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

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

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

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

  1230. 

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

  1232. 	 * the consumption of O2*/

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

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

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

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

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

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

  1239. 		aream=area;

  1240. 	}

  1241. 

  1242. 	//double maxOH,maxHO2;

  1243. 	//maxOH=0.0e0;

  1244. 	//maxHO2=0.0e0;

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

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

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

  1248. 	//	}

  1249. 	//}

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

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

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

  1253. 	//	}

  1254. 	//}

  1255. 

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

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

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

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

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

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

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

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

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

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

  1266. }

  1267. 

  1268. //

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

  1270. //{

  1271. //	double *ydata,*psidata;

  1272. //  	ydata    = N_VGetArrayPointer_OpenMP(y);

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

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

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

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

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

  1278. //		}

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

  1280. //	}

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

  1282. //}

  1283. //

  1284. //

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

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

  1287. ////  	double *ydata;

  1288. ////  	double *ydotdata;

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

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

  1291. ////

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

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

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

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

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

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

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

  1299. ////

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

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

  1302. ////	}

  1303. ////}