weiye
/
Droplet_LTORC
Beobachten 1
Favorisieren 0
Fork 0
Code Issues 0 Pull-Requests 0 Releases 0 Wiki Aktivität
Droplet Lagrangian Transient One-dimensional Reacting Code Implementation of both liquid and gas phase governing equations.
Du kannst nicht mehr als 25 Themen auswählen Themen müssen entweder mit einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.
1 Commit
1 Branch
1.7MB
Struktur: 504699f375
Branches Tags
${ item.name }
Erstelle Branch ${ searchTerm }
von „504699f375“
${ noResults }
Droplet_LTORC/residue.cpp

1423 Zeilen
42KB
Originalformat Blame Verlauf 

  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. //Locate species index:

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

  210. 	size_t index=0;

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

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

  213. 			index=k;

  214. 		}

  215. 	}

  216. 	return(index);

  217. }

  218. 

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

  220. 	double *iddata;

  221.   	N_VConst(ONE, *id);

  222.   	iddata   = N_VGetArrayPointer_OpenMP(*id);

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

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

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

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

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

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

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

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

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

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

  233. 		Rid(i)=ZERO;

  234. 		Pid(i)=ZERO;

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

  236. 		Mdotid(i)=ZERO;

  237. 	}	

  238. 	Mdotid(1)=ZERO;

  239. 	Rid(1)=ONE;

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

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

  242. 	if(data->constantPressure){

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

  244. 	}else{

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

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

  247. 	}

  248. 	if(data->dirichletInner){

  249. 		Tid(1)=ZERO; 

  250. 	}else{

  251. 		Tid(1)=ONE;

  252. 	}

  253. 	if(data->dirichletOuter){

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

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

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

  257. 			}

  258. 			Yid(1,k)=ZERO;

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

  260. 		}

  261. 	}else{

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

  263. 			Yid(1,k)=ZERO;

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

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

  266. 		}

  267. 	}

  268. 	return(0);

  269. }

  270. 

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

  272. 	switch (*i) {

  273. 		case 0:

  274. 			return(ONE);

  275. 		case 1:

  276. 			return(x);

  277. 		case 2:

  278. 			return(x*x);

  279. 		default:

  280. 			return(ONE);

  281. 	}

  282. }

  283. 

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

  285. 

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

  287. 

  288. 	size_t bufLen=10000;

  289. 	size_t nRows=0;

  290. 	size_t nColumns=nvar;

  291. 

  292. 	char buf[bufLen];

  293. 	char buf1[bufLen];

  294. 	char comment[1];

  295. 	char *ret;

  296. 

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

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

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

  300. 			nRows++;

  301. 		}

  302. 	}

  303. 	rewind(input);

  304. 

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

  306. 	double y[nRows*nColumns];

  307. 

  308. 	size_t i=0;

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

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

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

  312. 		}

  313. 		else{

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

  315. 			size_t j=0;

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

  317. 			j++;

  318. 			while(ret!=NULL){

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

  320. 				if(j<nColumns){

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

  322. 				}

  323. 				j++;

  324. 			}

  325. 			i++;

  326. 		}

  327. 	}

  328. 

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

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

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

  332. 		}

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

  334. 	}

  335. 	fclose(output);

  336. 

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

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

  339. 

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

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

  342. 	}

  343. 

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

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

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

  347. 	//}

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

  349. 	double dX0;

  350. 	int NN = nPts-2;

  351. 	

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

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

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

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

  356. 	}

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

  358. 

  359. 	gsl_interp_accel* acc;

  360. 	gsl_spline* spline;

  361. 	acc = gsl_interp_accel_alloc();

  362. 	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);

  363. 

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

  365. 

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

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

  368. 		}

  369. 

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

  371. 

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

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

  374. 		}

  375. 	}

  376. 	

  377. 	gsl_interp_accel_free(acc);

  378. 	gsl_spline_free(spline);

  379. 

  380. }

  381. 

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

  383. 	double mass=0.0e0;

  384. 	double rho;

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

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

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

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

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

  390. 	}

  391. 	return(mass);

  392. }

  393. 

  394. 

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

  396. 

  397. 	double rho,rhom;

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

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

  400. 	 * of kg. */

  401. 	psi(1)=ZERO;

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

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

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

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

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

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

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

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

  410. 	}

  411. 

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

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

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

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

  416. 	 * equations so that units match.*/

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

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

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

  420. 	}

  421. 	return(0);

  422. }

  423. 

  424. 

  425. int setInitialCondition(N_Vector* y, 

  426. 	    		N_Vector* ydot, 

  427.             		UserData data){

  428. 

  429. 	double* ydata;

  430. 	double* ydotdata;

  431. 	double* psidata;

  432. 	double* innerMassFractionsData, Rd, massDrop;

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

  434. 	double f=ZERO;

  435. 	double g=ZERO;

  436. 

  437. 	double perturb,rho;

  438. 	double epsilon=ZERO;

  439. 	int m,ier;

  440.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  441.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  442. 	innerMassFractionsData =  data->innerMassFractions;

  443. 	Rd = data->Rd;

  444. 	massDrop = data->massDrop;

  445. 

  446. 	//Mass of droplet

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

  448. 

  449. 	if(data->adaptiveGrid){

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

  451. 	}

  452. 	else{

  453.   		psidata  = data->uniformGrid;

  454. 	}

  455. 

  456. 	m=data->metric;

  457. 

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

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

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

  461. 	}

  462. 

  463. 	//Define Grid:

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

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

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

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

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

  469. 		}else{

  470. 			if(i==1){

  471. 				R(i)=ZERO;

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

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

  474. 			}else{

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

  476. 			}

  477. 		}

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

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

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

  481. 		}

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

  483. 	}	

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

  485. 

  486. 	double Tmax;

  487. 	double Tmin=data->initialTemperature;

  488. 	double w=data->mixingWidth;

  489. 	double YN2=ZERO;

  490. 	double YO2=ZERO;

  491. 	double YFuel,YOxidizer,sum;

  492. 

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

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

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

  496. 	//}

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

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

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

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

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

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

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

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

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

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

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

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

  509. 	//		}

  510. 	//	}

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

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

  513. 	//	}

  514. 	//}

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

  516. 	//	FILE* input;

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

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

  519. 	//		fclose(input);

  520. 	//	}

  521. 	//	else{

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

  523. 	//		return(-1);

  524. 	//	}

  525. 	//}

  526. 

  527. 	initializePsiGrid(ydata,psidata,data);

  528. 

  529. 	if(data->adaptiveGrid){

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

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

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

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

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

  535. 	//	}

  536. 	//	else{

  537. 	//	}

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

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

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

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

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

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

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

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

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

  547. 		}

  548. 	}

  549. 	else{

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

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

  552. 		double psiNew[data->npts];

  553. 		

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

  555. 		psiNew[0] = 0;

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

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

  558. 		}

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

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

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

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

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

  564. 

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

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

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

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

  569. 		//}

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

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

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

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

  574. 	}

  575. 

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

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

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

  579. 	}

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

  581. 		/*Premixed Combustion: Equilibrium products comprise ignition

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

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

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

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

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

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

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

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

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

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

  592. 			}

  593. 		}

  594. 		if(data->dirichletOuter){

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

  596. 		}

  597. 	}

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

  599. 		FILE* input;

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

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

  602. 			fclose(input);

  603. 		}

  604. 		else{

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

  606. 			return(-1);

  607. 		}

  608. 		initializePsiGrid(ydata,psidata,data);

  609. 	}

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

  611. 		FILE* input;

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

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

  614. 			fclose(input);

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

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

  617. 			return(0);

  618. 		}

  619. 		else{

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

  621. 			return(-1);

  622. 		}

  623. 	}

  624. 

  625. 	if(data->reflectProblem){

  626. 		double temp;

  627. 		int j=1;

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

  629. 			temp=T(j);

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

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

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

  633. 				temp=Y(j,k);

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

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

  636. 			}

  637. 			j=j+1;

  638. 		}

  639. 	}

  640. 

  641. 	/*Floor small values to zero*/

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

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

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

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

  646. 			}

  647. 		}

  648. 	}

  649. 

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

  651. 	if(data->adaptiveGrid){

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

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

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

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

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

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

  658. 	}

  659. 

  660. 	/*Ensure consistent boundary conditions*/	

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

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

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

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

  665. 	//}

  666. 

  667. 	return(0);

  668. }

  669. 

  670. inline double Qdot(double* t, 

  671. 	    	   double* x, 

  672. 	    	   double* ignTime, 

  673. 	    	   double* kernelSize, 

  674. 	    	   double* maxQdot){

  675. 	double qdot;

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

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

  678. 			qdot=(*maxQdot);

  679. 		}

  680. 		else{

  681. 			qdot=0.0e0;

  682. 		}

  683. 	}else{

  684. 		qdot=0.0e0;

  685. 	}

  686. 	return(qdot);

  687. }

  688. 

  689. inline void setGas(UserData data, 

  690. 	    double *ydata, 

  691. 	    size_t gridPoint){

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

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

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

  695. }

  696. 

  697. void getTransport(UserData data, 

  698. 		  double *ydata, 

  699. 		  size_t gridPoint,

  700. 		  double *rho,

  701. 		  double *lambda,

  702. 		  double YV[]){

  703. 

  704. 	double YAvg[data->nsp],

  705. 		 XLeft[data->nsp],

  706. 		 XRight[data->nsp],

  707. 		 gradX[data->nsp];

  708. 

  709. 	setGas(data,ydata,gridPoint);

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

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

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

  713.                                                                              

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

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

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

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

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

  719. 	}                                                                    

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

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

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

  723. 

  724. 

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

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

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

  728. 

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

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

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

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

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

  734. }

  735. 

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

  737. 

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

  739. 

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

  741.   	

  742. 	UserData data;

  743. 	data = (UserData)user_data;

  744. 	size_t npts=data->npts;

  745. 	size_t nsp=data->nsp;

  746. 	size_t k_bath = data->k_bath;

  747. 	size_t k_drop = data->k_drop;

  748. 

  749. 	ydata 	= N_VGetArrayPointer_OpenMP(y); 

  750. 	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 

  751. 	resdata = N_VGetArrayPointer_OpenMP(res);

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

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

  754. 	}else{

  755. 		psidata = data->uniformGrid;

  756. 	}

  757. 

  758. 	innerMassFractionsData =  data->innerMassFractions;

  759. 

  760. 	/* Grid stencil:*/

  761. 

  762. 	/*-------|---------*---------|---------*---------|-------*/

  763. 	/*-------|---------*---------|---------*---------|-------*/

  764. 	/*-------|---------*---------|---------*---------|-------*/

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

  766. 	/*-------|---------*---------|---------*---------|-------*/

  767. 	/*-------|---------*---------|---------*---------|-------*/

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

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

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

  771. 

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

  773. 	 * calculated values:

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

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

  776. 	 * m 		: exponent that determines geometry;

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

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

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

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

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

  782. 	 * Tav		: average temperature between two points.

  783. 	 * Pav		: average pressure between two points.

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

  785. 	 * Xgradhalf	: mole fraction gradient at j.

  786. 	 * Cpb		: mass based bulk specific heat.

  787. 	 * tranTerm	: transient terms.

  788. 	 * advTerm	: advection terms.

  789. 	 * diffTerm	: diffusion terms.

  790. 	 * srcTerm	: source terms.

  791. 	 */

  792. 

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

  794. 		 rho, 

  795. 		 rhophalf, lambdaphalf, YVphalf[nsp],

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

  797. 		 tranTerm, diffTerm,    srcTerm, advTerm,

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

  799. 	

  800. 	/*Aliases for difference coefficients:*/

  801. 	double cendfm, cendfc, cendfp;

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

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

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

  805. 	//double mass, mdotIn;

  806. 	double mass, massDrop;

  807. 	double sum, sum1, sum2, sum3;

  808. 

  809. 	size_t j,k;

  810. 	int m;

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

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

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

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

  815. 

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

  817. 	setGas(data,ydata,1);

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

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

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

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

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

  823. 	Cpl=  5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2)  //J/kg/K for liquid water

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

  825. 	

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

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

  828. 	dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01;	//J/kg	latent heat of vaporization of water

  829. 	rhol = 997.0;	

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

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

  832. 

  833. 	/*******************************************************************/

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

  835. 

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

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

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

  839. 	areamhalfsq= areamhalf*areamhalf;                                   

  840.                                                                              

  841. 	/*******************************************************************/

  842. 

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

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

  845. 

  846. 

  847. 	/*Mass:*/                                                            

  848. 	if(data->quasiSteady){

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

  850. 	}else{

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

  852. 	}

  853.                                                                         

  854. 	/*Species:*/                                                    

  855. 	sum=ZERO;

  856. 	

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

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

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

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

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

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

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

  864. 

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

  866. 

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

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

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

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

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

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

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

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

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

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

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

  878. 			//}

  879. 			//else{

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

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

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

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

  884. 			//	 * fix this?*/

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

  886. 			//}

  887. 		}

  888. 	}                                                                 

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

  890. 	

  891. 	/*Energy:*/                                                     

  892. 	if(data->dirichletInner){

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

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

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

  896. 	}else{

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

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

  899. 		//	/ (massDrop * Cpl);

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

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

  902. 			/ (massDrop * Cpl);

  903. 

  904. 

  905. 

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

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

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

  909. 	}

  910. 

  911. 	/*Pressure:*/                                                        

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

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

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

  915.                                                                              

  916. 		/*evaluate various mesh differences*///                      

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

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

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

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

  921. 		/***********************************///                      

  922.                                                                              

  923. 		/*evaluate various central difference coefficients*/       

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

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

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

  927. 		/**************************************************/

  928.                                                                    

  929.                                                                   

  930. 		/*evaluate properties at j*************************/

  931. 		setGas(data,ydata,j);

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

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

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

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

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

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

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

  939.                                                                            

  940. 		/*evaluate properties at p*************************/

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

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

  943. 		areaphalfsq= areaphalf*areaphalf;		    

  944. 		/**************************************************///     

  945. 

  946. 		/*Evaporating Mass*/

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

  948. 

  949. 		/*Mass:*/

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

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

  952. 

  953. 		/*Energy:*/

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

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

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

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

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

  959. 		/*Notes: 

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

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

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

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

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

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

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

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

  968. 		 * units.

  969. 		 * */

  970. 

  971. 		//enthalpy formulation:

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

  973. 		tranTerm = Tdot(j);

  974. 		sum=ZERO;

  975. 		sum1=ZERO;

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

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

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

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

  980. 		}

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

  982. 		sum1=sum1*Cantera::GasConstant;

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

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

  985. 			     /(dpsiav*Cpb) )

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

  987. 			              +cendfc*T(j)

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

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

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

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

  992. 			+advTerm 

  993. 			-diffTerm

  994. 			+srcTerm;

  995. 

  996. 	//	//energy formulation:

  997. 	//	tranTerm = Tdot(j);

  998. 	//	sum=ZERO;

  999. 	//	sum1=ZERO;

  1000. 	//	sum2=ZERO;

  1001. 	//	sum3=ZERO;

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

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

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

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

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

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

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

  1009. 	//	}

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

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

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

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

  1014. 	//		     /(dpsiav*Cvb) )

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

  1016. 	//		              +cendfc*T(j)

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

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

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

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

  1021. 	//	Tres(j)= tranTerm

  1022. 	//		+advTerm 

  1023. 	//		-diffTerm

  1024. 	//		+srcTerm;

  1025. 

  1026. 		/*Species:*/

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

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

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

  1030. 

  1031. 		sum=ZERO;

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

  1033. 			if(k!=k_bath){

  1034. 				tranTerm  = Ydot(j,k);

  1035. 				diffTerm  = (YVphalf(k)*areaphalf

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

  1037. 				srcTerm   = wdot(k)

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

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

  1040. 

  1041. 				Yres(j,k)= tranTerm 

  1042. 					  +advTerm

  1043. 					  +diffTerm

  1044. 					  -srcTerm;

  1045. 

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

  1047. 			}

  1048. 		}

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

  1050. 

  1051. 		/*Pressure:*/

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

  1053. 

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

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

  1056. 		areamhalf=areaphalf;                                

  1057. 		areamhalfsq=areaphalfsq;                            

  1058. 		aream=area;

  1059. 		rhom=rho;

  1060. 		rhomhalf=rhophalf;                                  

  1061. 		lambdamhalf=lambdaphalf;                            

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

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

  1064. 		}                                                   

  1065. 		/**************************************************/

  1066. 		                                                    

  1067. 	}                                                           

  1068. 	/*******************************************************************///

  1069. 

  1070. 

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

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

  1073. 	setGas(data,ydata,npts);

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

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

  1076. 

  1077. 	/*Mass:*/                                                          

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

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

  1080. 

  1081. 	/*Energy:*/                                             

  1082. 	if(data->dirichletOuter){

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

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

  1085. 	}

  1086. 	else{

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

  1088. 	}

  1089.                                                                 

  1090. 	/*Species:*/                                            

  1091. 	sum=ZERO;

  1092. 	if(data->dirichletOuter){

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

  1094. 			if(k!=k_bath){

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

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

  1097. 			}

  1098. 		}                                                       

  1099. 	}

  1100. 	else{

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

  1102. 			if(k!=k_bath){

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

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

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

  1106. 			}

  1107. 		}                                                       

  1108. 	}

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

  1110.                                                                 

  1111. 	

  1112. 	/*Pressure:*/                                                      

  1113. 	if(data->constantPressure){

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

  1115. 	}

  1116. 	else{

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

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

  1119. 	}

  1120. 

  1121. 	/*Evaporating Mass*/

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

  1123. 

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

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

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

  1127. 	//	//}

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

  1129. 	//}

  1130. 

  1131. 	return(0);

  1132. 

  1133. }

  1134. 

  1135. void printSpaceTimeHeader(UserData data)

  1136. {

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

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

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

  1140. 	}

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

  1142. 

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

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

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

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

  1147. 	}

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

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

  1150. }

  1151. 

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

  1153. {

  1154. 	double *ydata,*psidata;

  1155.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1156. 

  1157. 	if(data->adaptiveGrid){

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

  1159. 	}else{

  1160.   		psidata  = data->uniformGrid;

  1161. 	}

  1162. 

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

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

  1165. 		if(i==0){

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

  1167. 		}

  1168. 		else{

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

  1170. 		}

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

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

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

  1174. 			}else{

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

  1176. 			}

  1177. 		}

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

  1179. 	}

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

  1181. }

  1182. 

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

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

  1185. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1186. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1187. 	if(data->adaptiveGrid){

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

  1189. 	}else{

  1190.   		psidata  = data->uniformGrid;

  1191. 	}

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

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

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

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

  1196. }

  1197. 

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

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

  1200. 	double t;

  1201. 	if(data->adaptiveGrid){

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

  1203. 	}else{

  1204.   		psidata  = data->uniformGrid;

  1205. 	}

  1206. 	ydata = N_VGetArrayPointer_OpenMP(*y);

  1207. 	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

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

  1209. 	data->tNow=t;

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

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

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

  1213. 	if(data->adaptiveGrid){

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

  1215. 	}

  1216. }

  1217. 

  1218. void printGlobalHeader(UserData data)

  1219. {

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

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

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

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

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

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

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

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

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

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

  1230. }

  1231. //

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

  1233. //{

  1234. //	double *ydotdata,*psidata;

  1235. //  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);

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

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

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

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

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

  1241. //		}

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

  1243. //	}

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

  1245. //}

  1246. //

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

  1248. {

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

  1250. 	innerMassFractionsData =  data->innerMassFractions;

  1251. 

  1252. 	if(data->adaptiveGrid){

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

  1254. 	}else{

  1255.   		psidata  = data->uniformGrid;

  1256. 	}

  1257. 

  1258. 	double TAvg, RAvg, YAvg, psiAvg;

  1259.   	ydata    = N_VGetArrayPointer_OpenMP(*y);

  1260.   	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);

  1261. 	TAvg=data->isotherm;

  1262. 	double sum=ZERO;

  1263. 	double dpsim,area,aream,drdt;

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

  1265. 

  1266. 

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

  1268. 	size_t j=1;

  1269. 	size_t jj=1;

  1270. 	size_t jjj=1;

  1271. 	double wdot[data->nsp];

  1272. 	double wdotMax=0.0e0;

  1273. 	double advTerm=0.0e0;

  1274. 	psiAvg=0.0e0;

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

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

  1277. 			j=j+1;

  1278. 		}

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

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

  1281. 			jj=jj+1;

  1282. 		}

  1283. 

  1284. 		setGas(data,ydata,jj);

  1285. 		Tu=T(jj);

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

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

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

  1289. 		gamma=Cpb/Cvb;

  1290. 		

  1291. 	}

  1292. 	else{

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

  1294. 			j=j+1;

  1295. 		}

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

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

  1298. 			jj=jj+1;

  1299. 		}

  1300. 

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

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

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

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

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

  1306. 		gamma=Cpb/Cvb;

  1307. 	}

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

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

  1310. 	}

  1311. 	else{

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

  1313. 	}

  1314. 

  1315. 	////Experimental burning speed calculation:	

  1316. 	//int nMax=0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1335. 	////		}

  1336. 	////	}

  1337. 	////}

  1338. 	////else{

  1339. 	////	psiAvg=0.0e0;

  1340. 	////}

  1341. 

  1342. 

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

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

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

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

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

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

  1349. 

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

  1351. 	 * the consumption of O2*/

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

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

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

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

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

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

  1358. 		aream=area;

  1359. 	}

  1360. 

  1361. 	//double maxOH,maxHO2;

  1362. 	//maxOH=0.0e0;

  1363. 	//maxHO2=0.0e0;

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

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

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

  1367. 	//	}

  1368. 	//}

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

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

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

  1372. 	//	}

  1373. 	//}

  1374. 

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

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

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

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

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

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

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

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

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

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

  1385. }

  1386. 

  1387. //

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

  1389. //{

  1390. //	double *ydata,*psidata;

  1391. //  	ydata    = N_VGetArrayPointer_OpenMP(y);

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

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

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

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

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

  1397. //		}

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

  1399. //	}

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

  1401. //}

  1402. //

  1403. //

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

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

  1406. ////  	double *ydata;

  1407. ////  	double *ydotdata;

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

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

  1410. ////

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

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

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

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

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

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

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

  1418. ////

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

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

  1421. ////	}

  1422. ////}