weiye
/
HRR
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
Add the HRR related function and parameters into LTORC
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.
10 Commits
1 Branch
154KB
Branch: master
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'master'
${ noResults }
HRR/residue.cpp

1334 lines
38KB
Raw Permalink 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. 

  745. 	/*define the heat release rate related parameters*/

  746. 	//double Tsp=298.0; 

  747. 	double HRR = 0 ; 

  748. 	double Hf[nsp]; 

  749. 	//double heatRR[npts]; 

  750. 

  751. 	//double Hf = 0 ;

  752. 

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

  754. 	double mass, mdotIn;

  755. 	double sum, sum1, sum2, sum3;

  756. 

  757. 	size_t j,k;

  758. 	int m;

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

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

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

  762. 

  763. 	/*Initialize the HRR data*/ 

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

  765. 		HRRdata(j) = 0 ;

  766. 	}

  767. 

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

  769. 	setGas(data,ydata,1);

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

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

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

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

  774. 

  775. 	/*******************************************************************/

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

  777. 

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

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

  780. 	areamhalfsq= areamhalf*areamhalf;                                   

  781.                                                                              

  782. 	/*******************************************************************/

  783. 

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

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

  786. 

  787. 	/*Mass:*/                                                            

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

  789. 

  790. 	/*Energy:*/                                                     

  791. 

  792. 	if(data->dirichletInner){

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

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

  795. 	}

  796. 	else{

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

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

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

  800. 	}

  801.                                                                         

  802. 	/*Species:*/                                                    

  803. 	sum=ZERO;

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

  805. 		if(k!=k_bath){

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

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

  808. 					  Y(1,k)-

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

  810. 			}

  811. 			else{

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

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

  814. 				/*The below flux boundary condition makes the

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

  816. 				 * fix this?*/

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

  818. 			}

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

  820. 		}

  821. 	}                                                                 

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

  823. 

  824.                                                                              

  825. 	/*Pressure:*/                                                        

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

  827. 

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

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

  830.                                                                              

  831. 		/*evaluate various mesh differences*///                      

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

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

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

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

  836. 		/***********************************///                      

  837.                                                                              

  838. 		/*evaluate various central difference coefficients*/       

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

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

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

  842. 		/**************************************************/

  843.                                                                    

  844.                                                                   

  845. 		/*evaluate properties at j*************************/

  846. 		setGas(data,ydata,j);

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

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

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

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

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

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

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

  854.                                                                            

  855. 		/*evaluate properties at p*************************/

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

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

  858. 		areaphalfsq= areaphalf*areaphalf;		    

  859. 		/**************************************************///     

  860. 

  861. 		/*Mass:*/

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

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

  864. 

  865. 		/*Energy:*/

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

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

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

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

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

  871. 		/*Notes: 

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

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

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

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

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

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

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

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

  880. 		 * units.

  881. 		 * */

  882. 

  883. 		//enthalpy formulation:

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

  885. 		sum=ZERO;

  886. 		sum1=ZERO;

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

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

  889. 			sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*HALF*(YVmhalf(k)+YVphalf(k));

  890. 		}

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

  892. 		sum1=sum1*Cantera::GasConstant;

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

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

  895. 			     /(dpsiav*Cpb) )

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

  897. 			              +cendfc*T(j)

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

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

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

  901. 		Tres(j)= tranTerm

  902. 			+advTerm 

  903. 			-diffTerm

  904. 			+srcTerm;

  905. 		

  906. 		

  907. 		/*Calculate the Heat Release Rate */

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

  909. 			Hf(k) = data->gas->Hf298SS(k-1);

  910. 			HRR = - wdot(k) * Hf(k) ; 

  911. 			HRRdata(j) = HRR + HRRdata(j); 

  912. 		}

  913. 

  914. 	//	//energy formulation:

  915. 	//	tranTerm = Tdot(j);

  916. 	//	sum=ZERO;

  917. 	//	sum1=ZERO;

  918. 	//	sum2=ZERO;

  919. 	//	sum3=ZERO;

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

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

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

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

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

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

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

  927. 	//	}

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

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

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

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

  932. 	//		     /(dpsiav*Cvb) )

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

  934. 	//		              +cendfc*T(j)

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

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

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

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

  939. 	//	Tres(j)= tranTerm

  940. 	//		+advTerm 

  941. 	//		-diffTerm

  942. 	//		+srcTerm;

  943. 

  944. 		/*Species:*/

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

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

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

  948. 

  949. 		sum=ZERO;

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

  951. 			if(k!=k_bath){

  952. 				tranTerm  = Ydot(j,k);

  953. 				diffTerm  = (YVphalf(k)*areaphalf

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

  955. 				srcTerm   = wdot(k)

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

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

  958. 

  959. 				Yres(j,k)= tranTerm 

  960. 					  +advTerm

  961. 					  +diffTerm

  962. 					  -srcTerm;

  963. 

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

  965. 			}

  966. 		}

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

  968. 

  969. 		/*Pressure:*/

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

  971. 

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

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

  974. 		areamhalf=areaphalf;                                

  975. 		areamhalfsq=areaphalfsq;                            

  976. 		aream=area;

  977. 		rhom=rho;

  978. 		rhomhalf=rhophalf;                                  

  979. 		lambdamhalf=lambdaphalf;                            

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

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

  982. 		}                                                   

  983. 		/**************************************************/

  984. 		                                                    

  985. 	}                                                           

  986. 	/*******************************************************************///

  987. 

  988. 

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

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

  991. 	setGas(data,ydata,npts);

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

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

  994. 

  995. 	/*Mass:*/                                                          

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

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

  998. 

  999. 	/*Energy:*/                                             

  1000. 	if(data->dirichletOuter){

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

  1002. 	}

  1003. 	else{

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

  1005. 	}

  1006.                                                                 

  1007. 	/*Species:*/                                            

  1008. 	sum=ZERO;

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

  1010. 		if(k!=k_bath){

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

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

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

  1014. 		}

  1015. 	}                                                       

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

  1017.                                                                 

  1018. 	

  1019. 	/*Pressure:*/                                                      

  1020. 	if(data->constantPressure){

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

  1022. 	}

  1023. 	else{

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

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

  1026. 	}

  1027. 

  1028. 

  1029. 	/*******************************************************************/

  1030. 

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

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

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

  1034. 	//	//}

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

  1036. 	//}

  1037. 

  1038. 

  1039. 	return(0);

  1040. 

  1041. }

  1042. 

  1043. void printSpaceTimeHeader(UserData data)

  1044. {

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

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

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

  1048. 	}

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

  1050. 

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

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

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

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

  1055. 	}

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

  1057. 	fprintf((data->output), "%15s\t","HRR(J)");

  1058. 	fprintf((data->output), "%15s\n","rho(kg/m3)");

  1059. 

  1060. }

  1061. 

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

  1063. {

  1064. 	double *ydata,*psidata;

  1065. 	double rho;

  1066.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1067. 

  1068. 	if(data->adaptiveGrid){

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

  1070. 	}else{

  1071.   		psidata  = data->uniformGrid;

  1072. 	}

  1073. 

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

  1075. 		setGas(data,ydata,i+1) ; 

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

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

  1078. 		if(i==0){

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

  1080. 		}

  1081. 		else{

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

  1083. 		}

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

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

  1086. 		}

  1087. 		fprintf(output,"%15.6e\t",data->HRRdata[i]);

  1088. 		fprintf(output,"%15.9e\t",rho);

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

  1090. 	}

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

  1092. }

  1093. 

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

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

  1096. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1097. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1098. 	if(data->adaptiveGrid){

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

  1100. 	}else{

  1101.   		psidata  = data->uniformGrid;

  1102. 	}

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

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

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

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

  1107. }

  1108. 

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

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

  1111. 	double t;

  1112. 	if(data->adaptiveGrid){

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

  1114. 	}else{

  1115.   		psidata  = data->uniformGrid;

  1116. 	}

  1117. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1118. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

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

  1120. 	data->tNow=t;

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

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

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

  1124. 	if(data->adaptiveGrid){

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

  1126. 	}

  1127. }

  1128. 

  1129. void printGlobalHeader(UserData data)

  1130. {

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

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

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

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

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

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

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

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

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

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

  1141. }

  1142. //

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

  1144. //{

  1145. //	double *ydotdata,*psidata;

  1146. //  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);

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

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

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

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

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

  1152. //		}

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

  1154. //	}

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

  1156. //}

  1157. //

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

  1159. {

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

  1161. 	innerMassFractionsData =  data->innerMassFractions;

  1162. 

  1163. 	if(data->adaptiveGrid){

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

  1165. 	}else{

  1166.   		psidata  = data->uniformGrid;

  1167. 	}

  1168. 

  1169. 	double TAvg, RAvg, YAvg, psiAvg;

  1170.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1171.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1172. 	TAvg=data->isotherm;

  1173. 	double sum=ZERO;

  1174. 	double dpsim,area,aream,drdt;

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

  1176. 

  1177. 

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

  1179. 	size_t j=1;

  1180. 	size_t jj=1;

  1181. 	size_t jjj=1;

  1182. 	double wdot[data->nsp];

  1183. 	double wdotMax=0.0e0;

  1184. 	double advTerm=0.0e0;

  1185. 	psiAvg=0.0e0;

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

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

  1188. 			j=j+1;

  1189. 		}

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

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

  1192. 			jj=jj+1;

  1193. 		}

  1194. 

  1195. 		setGas(data,ydata,jj);

  1196. 		Tu=T(jj);

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

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

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

  1200. 		gamma=Cpb/Cvb;

  1201. 		

  1202. 	}

  1203. 	else{

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

  1205. 			j=j+1;

  1206. 		}

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

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

  1209. 			jj=jj+1;

  1210. 		}

  1211. 

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

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

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

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

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

  1217. 		gamma=Cpb/Cvb;

  1218. 	}

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

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

  1221. 	}

  1222. 	else{

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

  1224. 	}

  1225. 

  1226. 	////Experimental burning speed calculation:	

  1227. 	//int nMax=0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1246. 	////		}

  1247. 	////	}

  1248. 	////}

  1249. 	////else{

  1250. 	////	psiAvg=0.0e0;

  1251. 	////}

  1252. 

  1253. 

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

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

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

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

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

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

  1260. 

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

  1262. 	 * the consumption of O2*/

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

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

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

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

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

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

  1269. 		aream=area;

  1270. 	}

  1271. 

  1272. 	//double maxOH,maxHO2;

  1273. 	//maxOH=0.0e0;

  1274. 	//maxHO2=0.0e0;

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

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

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

  1278. 	//	}

  1279. 	//}

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

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

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

  1283. 	//	}

  1284. 	//}

  1285. 

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

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

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

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

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

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

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

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

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

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

  1296. }

  1297. 

  1298. //

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

  1300. //{

  1301. //	double *ydata,*psidata;

  1302. //  	ydata    = N_VGetArrayPointer_OpenMP(y);

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

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

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

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

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

  1308. //		}

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

  1310. //	}

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

  1312. //}

  1313. //

  1314. //

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

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

  1317. ////  	double *ydata;

  1318. ////  	double *ydotdata;

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

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

  1321. ////

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

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

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

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

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

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

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

  1329. ////

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

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

  1332. ////	}

  1333. ////}