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.

1329 lines
38KB

  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. double maxGradPosition(const double* y, const size_t nt,
  12. const size_t nvar, const double* x, size_t nPts){
  13. double maxGradT=0.0e0;
  14. double gradT=0.0e0;
  15. double pos=0.0e0;
  16. size_t j,jm;
  17. for (size_t i = 1; i <nPts; i++) {
  18. j=i*nvar+nt;
  19. jm=(i-1)*nvar+nt;
  20. gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
  21. if (gradT>=maxGradT) {
  22. maxGradT=gradT;
  23. pos=x[i];
  24. }
  25. }
  26. return(pos);
  27. }
  28. int maxGradIndex(const double* y, const size_t nt,
  29. const size_t nvar, const double* x, size_t nPts){
  30. double maxGradT=0.0e0;
  31. double gradT=0.0e0;
  32. int pos=0.0e0;
  33. size_t j,jm;
  34. for (size_t i = 1; i <nPts; i++) {
  35. j=i*nvar+nt;
  36. jm=(i-1)*nvar+nt;
  37. gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
  38. if (gradT>=maxGradT) {
  39. maxGradT=gradT;
  40. pos=i;
  41. }
  42. }
  43. return(pos);
  44. }
  45. double maxCurvPosition(const double* y, const size_t nt,
  46. const size_t nvar, const double* x, size_t nPts){
  47. double maxCurvT=0.0e0;
  48. double gradTp=0.0e0;
  49. double gradTm=0.0e0;
  50. double curvT=0.0e0;
  51. double dx=0.0e0;
  52. double pos=0.0e0;
  53. size_t j,jm,jp;
  54. for (size_t i = 1; i <nPts-1; i++) {
  55. j=i*nvar+nt;
  56. jm=(i-1)*nvar+nt;
  57. jp=(i+1)*nvar+nt;
  58. gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
  59. gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
  60. dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
  61. curvT=(gradTp-gradTm)/dx;
  62. if (curvT>=maxCurvT) {
  63. maxCurvT=curvT;
  64. pos=x[i];
  65. }
  66. }
  67. return(pos);
  68. }
  69. int maxCurvIndex(const double* y, const size_t nt,
  70. const size_t nvar, const double* x, size_t nPts){
  71. double maxCurvT=0.0e0;
  72. double gradTp=0.0e0;
  73. double gradTm=0.0e0;
  74. double curvT=0.0e0;
  75. double dx=0.0e0;
  76. int pos=0;
  77. size_t j,jm,jp;
  78. for (size_t i = 1; i <nPts-1; i++) {
  79. j=i*nvar+nt;
  80. jm=(i-1)*nvar+nt;
  81. jp=(i+1)*nvar+nt;
  82. gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
  83. gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
  84. dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
  85. curvT=(gradTp-gradTm)/dx;
  86. if (curvT>=maxCurvT) {
  87. maxCurvT=curvT;
  88. pos=i;
  89. }
  90. }
  91. return(pos);
  92. }
  93. double isothermPosition(const double* y, const double T, const size_t nt,
  94. const size_t nvar, const double* x, const size_t nPts){
  95. double pos=x[nPts-1];
  96. size_t j;
  97. for (size_t i = 1; i <nPts; i++) {
  98. j=i*nvar+nt;
  99. if (y[j]<=T) {
  100. pos=x[i];
  101. break;
  102. }
  103. }
  104. return(pos);
  105. }
  106. void updateSolution(double* y, double* ydot, const size_t nvar,
  107. const double xOld[],const double xNew[],const size_t nPts){
  108. double ytemp[nPts],ydottemp[nPts];
  109. gsl_interp_accel* acc;
  110. gsl_spline* spline;
  111. acc = gsl_interp_accel_alloc();
  112. spline = gsl_spline_alloc(gsl_interp_steffen, nPts);
  113. gsl_interp_accel* accdot;
  114. gsl_spline* splinedot;
  115. accdot = gsl_interp_accel_alloc();
  116. splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);
  117. for (size_t j = 0; j < nvar; j++) {
  118. for (size_t i = 0; i < nPts; i++) {
  119. ytemp[i]=y[j+i*nvar];
  120. ydottemp[i]=ydot[j+i*nvar];
  121. }
  122. gsl_spline_init(spline,xOld,ytemp,nPts);
  123. gsl_spline_init(splinedot,xOld,ydottemp,nPts);
  124. for (size_t i = 0; i < nPts; i++) {
  125. y[j+i*nvar]=gsl_spline_eval(spline,xNew[i],acc);
  126. ydot[j+i*nvar]=gsl_spline_eval(splinedot,xNew[i],accdot);
  127. }
  128. }
  129. //Exploring "fixing" boundary conditions:
  130. //for (size_t j = 1; j < nvar; j++) {
  131. // //printf("%15.6e\t%15.6e\n", y[j],y[j+nvar]);
  132. // y[j]=y[j+nvar];
  133. // //y[j+(nPts-1)*nvar]=y[j+(nPts-2)*nvar];
  134. // //ydot[j+nvar]=ydot[j];
  135. //}
  136. //y[0]=0.0e0;
  137. gsl_interp_accel_free(acc);
  138. gsl_spline_free(spline);
  139. gsl_interp_accel_free(accdot);
  140. gsl_spline_free(splinedot);
  141. }
  142. //Locate bath gas:
  143. size_t BathGasIndex(UserData data){
  144. size_t index=0;
  145. double max1;
  146. double max=data->gas->massFraction(0);
  147. for(size_t k=1;k<data->nsp;k++){
  148. max1=data->gas->massFraction(k);
  149. if(max1>=max){
  150. max=max1;
  151. index=k;
  152. }
  153. }
  154. return(index+1);
  155. }
  156. //Locate Oxidizer:
  157. size_t oxidizerIndex(UserData data){
  158. size_t index=0;
  159. for(size_t k=1;k<data->nsp;k++){
  160. if(data->gas->speciesName(k-1)=="O2"){
  161. index=k;
  162. }
  163. }
  164. return(index);
  165. }
  166. //Locate OH:
  167. size_t OHIndex(UserData data){
  168. size_t index=0;
  169. for(size_t k=1;k<data->nsp;k++){
  170. if(data->gas->speciesName(k-1)=="OH"){
  171. index=k;
  172. }
  173. }
  174. return(index);
  175. }
  176. //Locate HO2:
  177. size_t HO2Index(UserData data){
  178. size_t index=0;
  179. for(size_t k=1;k<data->nsp;k++){
  180. if(data->gas->speciesName(k-1)=="HO2"){
  181. index=k;
  182. }
  183. }
  184. return(index);
  185. }
  186. int setAlgebraicVariables(N_Vector* id, UserData data){
  187. double *iddata;
  188. N_VConst(ONE, *id);
  189. iddata = N_VGetArrayPointer_OpenMP(*id);
  190. data->k_bath=BathGasIndex(data);
  191. data->k_oxidizer=oxidizerIndex(data);
  192. data->k_OH=OHIndex(data);
  193. data->k_HO2=HO2Index(data);
  194. printf("Oxidizer index: %lu\n",data->k_oxidizer);
  195. printf("Bath gas index:%lu\n",data->k_bath);
  196. for (size_t i = 1; i <=data->npts; i++) {
  197. /*Algebraic variables: indicated by ZERO.*/
  198. Rid(i)=ZERO;
  199. Pid(i)=ZERO;
  200. Yid(i,data->k_bath)=ZERO;
  201. }
  202. Rid(1)=ONE;
  203. Tid(1)=ZERO;
  204. Tid(data->npts)=ZERO;
  205. if(data->constantPressure){
  206. Pid(data->npts)=ONE;
  207. }
  208. else{
  209. Pid(data->npts)=ZERO;
  210. Rid(data->npts)=ONE;
  211. }
  212. for (size_t k = 1; k <=data->nsp; k++) {
  213. Yid(1,k)=ZERO;
  214. Yid(data->npts,k)=ZERO;
  215. }
  216. return(0);
  217. }
  218. inline double calc_area(double x,int* i){
  219. switch (*i) {
  220. case 0:
  221. return(ONE);
  222. case 1:
  223. return(x);
  224. case 2:
  225. return(x*x);
  226. default:
  227. return(ONE);
  228. }
  229. }
  230. void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts){
  231. FILE* output;output=fopen("test.dat","w");
  232. size_t bufLen=10000;
  233. size_t nRows=0;
  234. size_t nColumns=nvar;
  235. char buf[bufLen];
  236. char buf1[bufLen];
  237. char comment[1];
  238. char *ret;
  239. while (fgets(buf,bufLen, input)!=NULL){
  240. comment[0]=buf[0];
  241. if(strncmp(comment,"#",1)!=0){
  242. nRows++;
  243. }
  244. }
  245. rewind(input);
  246. printf("nRows: %ld\n", nRows);
  247. double y[nRows*nColumns];
  248. size_t i=0;
  249. while (fgets(buf,bufLen, input)!=NULL){
  250. comment[0]=buf[0];
  251. if(strncmp(comment,"#",1)==0){
  252. }
  253. else{
  254. ret=strtok(buf,"\t");
  255. size_t j=0;
  256. y[i*nColumns+j]=(double)(atof(ret));
  257. j++;
  258. while(ret!=NULL){
  259. ret=strtok(NULL,"\t");
  260. if(j<nColumns){
  261. y[i*nColumns+j]=(double)(atof(ret));
  262. }
  263. j++;
  264. }
  265. i++;
  266. }
  267. }
  268. for (i = 0; i < nRows; i++) {
  269. for (size_t j = 0; j < nColumns; j++) {
  270. fprintf(output, "%15.6e\t",y[i*nColumns+j]);
  271. }
  272. fprintf(output, "\n");
  273. }
  274. fclose(output);
  275. double xOld[nRows],xNew[nPts],ytemp[nPts];
  276. for (size_t j = 0; j < nRows; j++) {
  277. xOld[j]=y[j*nColumns+nr];
  278. }
  279. double dx=xOld[nRows-1]/((double)(nPts)-1.0e0);
  280. for (size_t j = 0; j < nPts; j++) {
  281. xNew[j]=(double)j*dx;
  282. }
  283. gsl_interp_accel* acc;
  284. gsl_spline* spline;
  285. acc = gsl_interp_accel_alloc();
  286. spline = gsl_spline_alloc(gsl_interp_steffen, nRows);
  287. for (size_t j = 0; j < nColumns; j++) {
  288. for (size_t k = 0; k < nRows; k++) {
  289. ytemp[k]=y[j+k*nColumns];
  290. }
  291. gsl_spline_init(spline,xOld,ytemp,nRows);
  292. for (size_t k = 1; k < nPts-1; k++) {
  293. ydata[j+k*nColumns]=gsl_spline_eval(spline,xNew[k],acc);
  294. }
  295. }
  296. gsl_interp_accel_free(acc);
  297. gsl_spline_free(spline);
  298. }
  299. double systemMass(double* ydata,UserData data){
  300. double mass=0.0e0;
  301. double rho;
  302. for (size_t i = 2; i <=data->npts; i++) {
  303. data->gas->setState_TPY(T(i), P(i), &Y(i,1));
  304. rho=data->gas->density();
  305. //psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&m);
  306. mass+=rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
  307. }
  308. return(mass);
  309. }
  310. int initializePsiGrid(double* ydata, double* psidata, UserData data){
  311. double rho;
  312. /*Create a psi grid that corresponds CONSISTENTLY to the spatial grid
  313. * "R" created above. Note that the Lagrangian variable psi has units
  314. * of kg. */
  315. psi(1)=ZERO;
  316. for (size_t i = 2; i <=data->npts; i++) {
  317. data->gas->setState_TPY(T(i), P(i), &Y(i,1));
  318. rho=data->gas->density();
  319. //psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&data->metric);
  320. psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
  321. }
  322. /*The mass of the entire system is the value of psi at the last grid
  323. * point. Normalize psi by this mass so that it varies from zero to
  324. * one. This makes psi dimensionless. So the mass needs to be
  325. * multiplied back in the approporiate places in the governing
  326. * equations so that units match.*/
  327. data->mass=psi(data->npts);
  328. for (size_t i = 1; i <=data->npts; i++) {
  329. psi(i)=psi(i)/data->mass;
  330. }
  331. return(0);
  332. }
  333. int setInitialCondition(N_Vector* y,
  334. N_Vector* ydot,
  335. UserData data){
  336. double* ydata;
  337. double* ydotdata;
  338. double* psidata;
  339. double* innerMassFractionsData;
  340. double f=ZERO;
  341. double g=ZERO;
  342. double perturb,rho;
  343. double epsilon=ZERO;
  344. int m,ier;
  345. ydata = N_VGetArrayPointer_OpenMP(*y);
  346. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  347. innerMassFractionsData = data->innerMassFractions;
  348. if(data->adaptiveGrid){
  349. psidata = data->grid->xOld;
  350. }
  351. else{
  352. psidata = data->uniformGrid;
  353. }
  354. m=data->metric;
  355. data->innerTemperature=data->initialTemperature;
  356. for (size_t k = 1; k <=data->nsp; k++) {
  357. innerMassFractionsData[k-1]=data->gas->massFraction(k-1);
  358. }
  359. //Define Grid:
  360. //double Rmin=1e-03*data->domainLength;
  361. double Rmin=0.0e0;
  362. double dR=(data->domainLength-Rmin)/((double)(data->npts)-1.0e0);
  363. double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0);
  364. for (size_t i = 1; i <=data->npts; i++) {
  365. if(data->metric==0){
  366. R(i)=Rmin+(double)((i-1)*dR);
  367. }else{
  368. if(i==1){
  369. R(i)=ZERO;
  370. }else if(i==2){
  371. R(i)=data->firstRadius*data->domainLength;
  372. }else{
  373. R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));
  374. }
  375. }
  376. T(i)=data->initialTemperature;
  377. for (size_t k = 1; k <=data->nsp; k++) {
  378. Y(i,k)=data->gas->massFraction(k-1); //Indexing different in Cantera
  379. }
  380. P(i)=data->initialPressure*Cantera::OneAtm;
  381. }
  382. R(data->npts)=data->domainLength;
  383. double Tmax;
  384. double Tmin=data->initialTemperature;
  385. double w=data->mixingWidth;
  386. double YN2=ZERO;
  387. double YO2=ZERO;
  388. double YFuel,YOxidizer,sum;
  389. //if(data->problemType==0){
  390. // data->gas->equilibrate("HP");
  391. // data->maxTemperature=data->gas->temperature();
  392. //}
  393. //else if(data->problemType==1){
  394. // /*Premixed Combustion: Equilibrium products comprise ignition
  395. // * kernel at t=0. The width of the kernel is "mixingWidth"
  396. // * shifted by "shift" from the center.*/
  397. // data->gas->equilibrate("HP");
  398. // Tmax=data->gas->temperature();
  399. // for (size_t i = 1; i <=data->npts; i++) {
  400. // g=HALF*(tanh((R(i)-data->shift)/w)+ONE); //increasing function of x
  401. // f=ONE-g; //decreasing function of x
  402. // T(i)=(Tmax-Tmin)*f+Tmin;
  403. // for (size_t k = 1; k <=data->nsp; k++) {
  404. // Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
  405. // }
  406. // }
  407. // if(data->dirichletOuter){
  408. // T(data->npts)=data->wallTemperature;
  409. // }
  410. //}
  411. //else if(data->problemType==2){
  412. // FILE* input;
  413. // if(input=fopen("initialCondition.dat","r")){
  414. // readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);
  415. // fclose(input);
  416. // }
  417. // else{
  418. // printf("file initialCondition.dat not found!\n");
  419. // return(-1);
  420. // }
  421. //}
  422. initializePsiGrid(ydata,psidata,data);
  423. if(data->adaptiveGrid){
  424. // if(data->problemType!=0){
  425. // data->grid->position=maxGradPosition(ydata, data->nt, data->nvar,
  426. // data->grid->xOld, data->npts);
  427. // //data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
  428. // // data->grid->xOld, data->npts);
  429. // }
  430. // else{
  431. // }
  432. if(data->problemType!=3){
  433. data->grid->position=0.0e0;
  434. double x=data->grid->position+data->gridOffset*data->grid->leastMove;
  435. printf("New grid center:%15.6e\n",x);
  436. ier=reGrid(data->grid, x);
  437. if(ier==-1)return(-1);
  438. updateSolution(ydata, ydotdata, data->nvar,
  439. data->grid->xOld,data->grid->x,data->npts);
  440. storeGrid(data->grid->x,data->grid->xOld,data->npts);
  441. }
  442. }
  443. else{
  444. double psiNew[data->npts];
  445. double dpsi=1.0e0/((double)(data->npts)-1.0e0);
  446. for (size_t i = 0; i < data->npts; i++) {
  447. psiNew[i]=(double)(i)*dpsi;
  448. }
  449. printf("Last point:%15.6e\n",psiNew[data->npts-1]);
  450. updateSolution(ydata, ydotdata, data->nvar,
  451. data->uniformGrid,psiNew,data->npts);
  452. storeGrid(psiNew,data->uniformGrid,data->npts);
  453. }
  454. if(data->problemType==0){
  455. data->gas->equilibrate("HP");
  456. data->maxTemperature=data->gas->temperature();
  457. }
  458. else if(data->problemType==1){
  459. /*Premixed Combustion: Equilibrium products comprise ignition
  460. * kernel at t=0. The width of the kernel is "mixingWidth"
  461. * shifted by "shift" from the center.*/
  462. data->gas->equilibrate("HP");
  463. Tmax=data->gas->temperature();
  464. for (size_t i = 1; i <=data->npts; i++) {
  465. g=HALF*(tanh((R(i)-data->shift)/w)+ONE); //increasing function of x
  466. f=ONE-g; //decreasing function of x
  467. T(i)=(Tmax-Tmin)*f+Tmin;
  468. for (size_t k = 1; k <=data->nsp; k++) {
  469. Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
  470. }
  471. }
  472. if(data->dirichletOuter){
  473. T(data->npts)=data->wallTemperature;
  474. }
  475. }
  476. else if(data->problemType==2){
  477. FILE* input;
  478. if(input=fopen("initialCondition.dat","r")){
  479. readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);
  480. fclose(input);
  481. }
  482. else{
  483. printf("file initialCondition.dat not found!\n");
  484. return(-1);
  485. }
  486. }
  487. else if(data->problemType==3){
  488. FILE* input;
  489. if(input=fopen("restart.bin","r")){
  490. readRestart(y, ydot, input, data);
  491. fclose(input);
  492. printf("Restart solution loaded!\n");
  493. printf("Problem starting at t=%15.6e\n",data->tNow);
  494. return(0);
  495. }
  496. else{
  497. printf("file restart.bin not found!\n");
  498. return(-1);
  499. }
  500. }
  501. if(data->reflectProblem){
  502. double temp;
  503. int j=1;
  504. while (data->npts+1-2*j>=0) {
  505. temp=T(j);
  506. T(j)=T(data->npts+1-j);
  507. T(data->npts+1-j)=temp;
  508. for (size_t k = 1; k <=data->nsp; k++) {
  509. temp=Y(j,k);
  510. Y(j,k)=Y(data->npts+1-j,k);
  511. Y(data->npts+1-j,k)=temp;
  512. }
  513. j=j+1;
  514. }
  515. }
  516. /*Floor small values to zero*/
  517. for (size_t i = 1; i <=data->npts; i++) {
  518. for (size_t k = 1; k <=data->nsp; k++) {
  519. if(fabs(Y(i,k))<=data->massFractionTolerance){
  520. Y(i,k)=0.0e0;
  521. }
  522. }
  523. }
  524. //Set grid to location of maximum curvature:
  525. if(data->adaptiveGrid){
  526. data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
  527. data->grid->x, data->npts);
  528. ier=reGrid(data->grid, data->grid->position);
  529. updateSolution(ydata, ydotdata, data->nvar,
  530. data->grid->xOld,data->grid->x,data->npts);
  531. storeGrid(data->grid->x,data->grid->xOld,data->npts);
  532. }
  533. /*Ensure consistent boundary conditions*/
  534. T(1)=T(2);
  535. for (size_t k = 1; k <=data->nsp; k++) {
  536. Y(1,k)=Y(2,k);
  537. Y(data->npts,k)=Y(data->npts-1,k);
  538. }
  539. return(0);
  540. }
  541. inline double Qdot(double* t,
  542. double* x,
  543. double* ignTime,
  544. double* kernelSize,
  545. double* maxQdot){
  546. double qdot;
  547. if(*x<=*kernelSize){
  548. if((*t)<=(*ignTime)){
  549. qdot=(*maxQdot);
  550. }
  551. else{
  552. qdot=0.0e0;
  553. }
  554. }else{
  555. qdot=0.0e0;
  556. }
  557. return(qdot);
  558. }
  559. inline void setGas(UserData data,
  560. double *ydata,
  561. size_t gridPoint){
  562. data->gas->setTemperature(T(gridPoint));
  563. data->gas->setMassFractions_NoNorm(&Y(gridPoint,1));
  564. data->gas->setPressure(P(gridPoint));
  565. }
  566. void getTransport(UserData data,
  567. double *ydata,
  568. size_t gridPoint,
  569. double *rho,
  570. double *lambda,
  571. double YV[]){
  572. double YAvg[data->nsp],
  573. XLeft[data->nsp],
  574. XRight[data->nsp],
  575. gradX[data->nsp];
  576. setGas(data,ydata,gridPoint);
  577. data->gas->getMoleFractions(XLeft);
  578. setGas(data,ydata,gridPoint+1);
  579. data->gas->getMoleFractions(XRight);
  580. for (size_t k = 1; k <=data->nsp; k++) {
  581. YAvg(k)=HALF*(Y(gridPoint,k)+
  582. Y(gridPoint+1,k));
  583. gradX(k)=(XRight(k)-XLeft(k))/
  584. (R(gridPoint+1)-R(gridPoint));
  585. }
  586. double TAvg = HALF*(T(gridPoint)+T(gridPoint+1));
  587. double gradT=(T(gridPoint+1)-T(gridPoint))/
  588. (R(gridPoint+1)-R(gridPoint));
  589. data->gas->setTemperature(TAvg);
  590. data->gas->setMassFractions_NoNorm(YAvg);
  591. data->gas->setPressure(P(gridPoint));
  592. *rho=data->gas->density();
  593. *lambda=data->trmix->thermalConductivity();
  594. data->trmix->getSpeciesFluxes(1,&gradT,data->nsp,
  595. gradX,data->nsp,YV);
  596. //setGas(data,ydata,gridPoint);
  597. }
  598. int residue(double t,
  599. N_Vector y,
  600. N_Vector ydot,
  601. N_Vector res,
  602. void *user_data){
  603. /*Declare and fetch nvectors and user data:*/
  604. double *ydata, *ydotdata, *resdata, *psidata, *innerMassFractionsData;
  605. UserData data;
  606. data = (UserData)user_data;
  607. size_t npts=data->npts;
  608. size_t nsp=data->nsp;
  609. size_t k_bath = data->k_bath;
  610. ydata = N_VGetArrayPointer_OpenMP(y);
  611. ydotdata= N_VGetArrayPointer_OpenMP(ydot);
  612. resdata = N_VGetArrayPointer_OpenMP(res);
  613. if(data->adaptiveGrid==1){
  614. psidata = data->grid->x;
  615. }else{
  616. psidata = data->uniformGrid;
  617. }
  618. innerMassFractionsData = data->innerMassFractions;
  619. /* Grid stencil:*/
  620. /*-------|---------*---------|---------*---------|-------*/
  621. /*-------|---------*---------|---------*---------|-------*/
  622. /*-------|---------*---------|---------*---------|-------*/
  623. /*-------m-------mhalf-------j-------phalf-------p-------*/
  624. /*-------|---------*---------|---------*---------|-------*/
  625. /*-------|---------*---------|---------*---------|-------*/
  626. /*-------|<=======dxm=======>|<=======dxp=======>|-------*/
  627. /*-------|---------*<======dxav=======>*---------|-------*/
  628. /*-------|<================dxpm=================>|-------*/
  629. /* Various variables defined for book-keeping and storing previously
  630. * calculated values:
  631. * rho : densities at points m, mhalf, j, p, and phalf.
  632. * area : the matric at points m, mhalf, j, p, and phalf.
  633. * m : exponent that determines geometry;
  634. * lambda : thermal conductivities at mhalf and phalf.
  635. * mdot : mass flow rate at m, j, and p.
  636. * X : mole fractions at j and p.
  637. * YV : diffusion fluxes at mhalf and phalf.
  638. * Tgrad : temperature gradient at mhalf and phalf.
  639. * Tav : average temperature between two points.
  640. * Pav : average pressure between two points.
  641. * Yav : average mass fractions between two points.
  642. * Xgradhalf : mole fraction gradient at j.
  643. * Cpb : mass based bulk specific heat.
  644. * tranTerm : transient terms.
  645. * advTerm : advection terms.
  646. * diffTerm : diffusion terms.
  647. * srcTerm : source terms.
  648. */
  649. double rhomhalf, rhom, lambdamhalf, YVmhalf[nsp],
  650. rho,
  651. rhophalf, lambdaphalf, YVphalf[nsp],
  652. Cpb, Cvb, Cp[nsp], wdot[nsp], enthalpy[nsp], energy[nsp],
  653. tranTerm, diffTerm, srcTerm, advTerm,
  654. area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq;
  655. /*Aliases for difference coefficients:*/
  656. double cendfm, cendfc, cendfp;
  657. /*Aliases for various grid spacings:*/
  658. double dpsip, dpsiav, dpsipm, dpsim, dpsimm;
  659. /*define the heat release rate related parameters*/
  660. //double Tsp=298.0;
  661. double HRR = 0 ;
  662. double Hf[nsp];
  663. //double heatRR[npts];
  664. //double Hf = 0 ;
  665. dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;
  666. double mass, mdotIn;
  667. double sum, sum1, sum2, sum3;
  668. size_t j,k;
  669. int m;
  670. m=data->metric; //Unitless
  671. mass=data->mass; //Units: kg
  672. mdotIn=data->mdot*calc_area(R(npts),&m); //Units: kg/s
  673. /*Initialize the HRR data*/
  674. for (j=1; j<= npts ; j++) {
  675. HRRdata(j) = 0 ;
  676. }
  677. // /*evaluate properties at j=1*************************/
  678. setGas(data,ydata,1);
  679. rhom=data->gas->density();
  680. Cpb=data->gas->cp_mass(); //J/kg/K
  681. Cvb=data->gas->cv_mass(); //J/kg/K
  682. aream= calc_area(R(1),&m);
  683. /*******************************************************************/
  684. /*Calculate values at j=2's m and mhalf*****************************/
  685. getTransport(data, ydata, 1, &rhomhalf,&lambdamhalf,YVmhalf);
  686. areamhalf= calc_area(HALF*(R(1)+R(2)),&m);
  687. areamhalfsq= areamhalf*areamhalf;
  688. /*******************************************************************/
  689. /*Fill up res with left side (center) boundary conditions:**********/
  690. /*We impose zero fluxes at the center:*/
  691. /*Mass:*/
  692. Rres(1)=Rdot(1);
  693. /*Energy:*/
  694. if(data->dirichletInner){
  695. //Tres(1)=Tdot(1);
  696. Tres(1)=T(1)-data->innerTemperature;
  697. }
  698. else{
  699. Tres(1)=T(2)-T(1);
  700. //Tres(1)=Tdot(1) - (Pdot(1)/(rhom*Cpb))
  701. // +(double)(data->metric+1)*(rhomhalf*lambdamhalf*areamhalfsq*(T(2)-T(1))/psi(2)-psi(1));
  702. }
  703. /*Species:*/
  704. sum=ZERO;
  705. for (k = 1; k <=nsp; k++) {
  706. if(k!=k_bath){
  707. if(fabs(mdotIn)>1e-14){
  708. Yres(1,k)=innerMassFractionsData[k-1]-
  709. Y(1,k)-
  710. (YVmhalf(k)*areamhalf)/mdotIn;
  711. }
  712. else{
  713. //Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];
  714. Yres(1,k)=Y(2,k)-Y(1,k);
  715. /*The below flux boundary condition makes the
  716. * problem more prone to diverge. How does one
  717. * fix this?*/
  718. //Yres(1,k)=YVmhalf(k);
  719. }
  720. sum=sum+Y(1,k);
  721. }
  722. }
  723. Yres(1,k_bath)=ONE-sum-Y(1,k_bath);
  724. /*Pressure:*/
  725. Pres(1)=P(2)-P(1);
  726. /*Fill up res with governing equations at inner points:*************/
  727. for (j = 2; j < npts; j++) {
  728. /*evaluate various mesh differences*///
  729. dpsip = (psi(j+1) - psi(j) )*mass;
  730. dpsim = (psi(j) - psi(j-1))*mass;
  731. dpsiav = HALF*(psi(j+1) - psi(j-1))*mass;
  732. dpsipm = (psi(j+1) - psi(j-1))*mass;
  733. /***********************************///
  734. /*evaluate various central difference coefficients*/
  735. cendfm = - dpsip / (dpsim*dpsipm);
  736. cendfc = (dpsip-dpsim) / (dpsip*dpsim);
  737. cendfp = dpsim / (dpsip*dpsipm);
  738. /**************************************************/
  739. /*evaluate properties at j*************************/
  740. setGas(data,ydata,j);
  741. rho=data->gas->density(); //kg/m^3
  742. Cpb=data->gas->cp_mass(); //J/kg/K
  743. Cvb=data->gas->cv_mass(); //J/kg/K
  744. data->gas->getNetProductionRates(wdot); //kmol/m^3
  745. data->gas->getEnthalpy_RT(enthalpy); //unitless
  746. data->gas->getCp_R(Cp); //unitless
  747. area = calc_area(R(j),&m); //m^2
  748. /*evaluate properties at p*************************/
  749. getTransport(data, ydata, j, &rhophalf,&lambdaphalf,YVphalf);
  750. areaphalf= calc_area(HALF*(R(j)+R(j+1)),&m);
  751. areaphalfsq= areaphalf*areaphalf;
  752. /**************************************************///
  753. /*Mass:*/
  754. /* ∂r/∂ψ = 1/ρA */
  755. Rres(j)=((R(j)-R(j-1))/dpsim)-(TWO/(rhom*aream+rho*area));
  756. /*Energy:*/
  757. /* ∂T/∂t = - ṁ(∂T/∂ψ)
  758. * + (1/cₚ)(∂/∂ψ)(λρA²∂T/∂ψ)
  759. * - (A/cₚ) ∑ YᵢVᵢcₚᵢ(∂T/∂ψ)
  760. * - (1/ρcₚ)∑ ώᵢhᵢ
  761. * + (1/ρcₚ)(∂P/∂t) */
  762. /*Notes:
  763. * λ has units J/m/s/K.
  764. * YᵢVᵢ has units kg/m^2/s.
  765. * hᵢ has units J/kmol, so we must multiply the enthalpy
  766. * defined above (getEnthalpy_RT) by T (K) and the gas constant
  767. * (J/kmol/K) to get the right units.
  768. * cₚᵢ has units J/kg/K, so we must multiply the specific heat
  769. * defined above (getCp_R) by the gas constant (J/kmol/K) and
  770. * divide by the molecular weight (kg/kmol) to get the right
  771. * units.
  772. * */
  773. //enthalpy formulation:
  774. tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));
  775. sum=ZERO;
  776. sum1=ZERO;
  777. for (k = 1; k <=nsp; k++) {
  778. sum=sum+wdot(k)*enthalpy(k);
  779. sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*HALF*(YVmhalf(k)+YVphalf(k));
  780. }
  781. sum=sum*Cantera::GasConstant*T(j);
  782. sum1=sum1*Cantera::GasConstant;
  783. diffTerm =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
  784. -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
  785. /(dpsiav*Cpb) )
  786. -(sum1*area*(cendfp*T(j+1)
  787. +cendfc*T(j)
  788. +cendfm*T(j-1))/Cpb);
  789. srcTerm = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cpb);
  790. advTerm = (mdotIn*(T(j)-T(j-1))/dpsim);
  791. Tres(j)= tranTerm
  792. +advTerm
  793. -diffTerm
  794. +srcTerm;
  795. /*Calculate the Heat Release Rate */
  796. for(size_t k = 1; k <= nsp; k++) {
  797. Hf(k) = data->gas->Hf298SS(k-1);
  798. HRR = - wdot(k) * Hf(k) ;
  799. HRRdata(j) = HRR + HRRdata(j);
  800. }
  801. // //energy formulation:
  802. // tranTerm = Tdot(j);
  803. // sum=ZERO;
  804. // sum1=ZERO;
  805. // sum2=ZERO;
  806. // sum3=ZERO;
  807. // for (k = 1; k <=nsp; k++) {
  808. // energy(k)=enthalpy(k)-ONE;
  809. // sum=sum+wdot(k)*energy(k);
  810. // sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
  811. // *HALF*(YVmhalf(k)+YVphalf(k));
  812. // sum2=sum2+(YVmhalf(k)/data->gas->molecularWeight(k-1));
  813. // sum3=sum3+(YVphalf(k)/data->gas->molecularWeight(k-1));
  814. // }
  815. // sum=sum*Cantera::GasConstant*T(j);
  816. // sum1=sum1*Cantera::GasConstant;
  817. // diffTerm =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
  818. // -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
  819. // /(dpsiav*Cvb) )
  820. // -(sum1*area*(cendfp*T(j+1)
  821. // +cendfc*T(j)
  822. // +cendfm*T(j-1))/Cvb);
  823. // srcTerm = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cvb);
  824. // advTerm = (mdotIn*(T(j)-T(j-1))/dpsim);
  825. // advTerm = advTerm + (Cantera::GasConstant*T(j)*area/Cvb)*((sum3-sum2)/dpsiav);
  826. // Tres(j)= tranTerm
  827. // +advTerm
  828. // -diffTerm
  829. // +srcTerm;
  830. /*Species:*/
  831. /* ∂Yᵢ/∂t = - ṁ(∂Yᵢ/∂ψ)
  832. * - (∂/∂ψ)(AYᵢVᵢ)
  833. * + (ώᵢWᵢ/ρ) */
  834. sum=ZERO;
  835. for (k = 1; k <=nsp; k++) {
  836. if(k!=k_bath){
  837. tranTerm = Ydot(j,k);
  838. diffTerm = (YVphalf(k)*areaphalf
  839. -YVmhalf(k)*areamhalf)/dpsiav;
  840. srcTerm = wdot(k)
  841. *(data->gas->molecularWeight(k-1))/rho;
  842. advTerm = (mdotIn*(Y(j,k)-Y(j-1,k))/dpsim);
  843. Yres(j,k)= tranTerm
  844. +advTerm
  845. +diffTerm
  846. -srcTerm;
  847. sum=sum+Y(j,k);
  848. }
  849. }
  850. Yres(j,k_bath)=ONE-sum-Y(j,k_bath);
  851. /*Pressure:*/
  852. Pres(j) = P(j+1)-P(j);
  853. /*Assign values evaluated at p and phalf to m
  854. * and mhalf to save some cpu cost:****************/
  855. areamhalf=areaphalf;
  856. areamhalfsq=areaphalfsq;
  857. aream=area;
  858. rhom=rho;
  859. rhomhalf=rhophalf;
  860. lambdamhalf=lambdaphalf;
  861. for (k = 1; k <=nsp; k++) {
  862. YVmhalf(k)=YVphalf(k);
  863. }
  864. /**************************************************/
  865. }
  866. /*******************************************************************///
  867. /*Fill up res with right side (wall) boundary conditions:***********/
  868. /*We impose zero fluxes at the wall:*/
  869. setGas(data,ydata,npts);
  870. rho=data->gas->density();
  871. area = calc_area(R(npts),&m);
  872. /*Mass:*/
  873. dpsim=(psi(npts)-psi(npts-1))*mass;
  874. Rres(npts)=((R(npts)-R(npts-1))/dpsim)-(TWO/(rhom*aream+rho*area));
  875. /*Energy:*/
  876. if(data->dirichletOuter){
  877. Tres(npts)=T(npts)-data->wallTemperature;
  878. }
  879. else{
  880. Tres(npts)=T(npts)-T(npts-1);
  881. }
  882. /*Species:*/
  883. sum=ZERO;
  884. for (k = 1; k <=nsp; k++) {
  885. if(k!=k_bath){
  886. //Yres(npts,k)=Y(npts,k)-Y(npts-1,k);
  887. Yres(npts,k)=YVmhalf(k);
  888. sum=sum+Y(npts,k);
  889. }
  890. }
  891. Yres(npts,k_bath)=ONE-sum-Y(npts,k_bath);
  892. /*Pressure:*/
  893. if(data->constantPressure){
  894. Pres(npts)=Pdot(npts)-data->dPdt;
  895. }
  896. else{
  897. Pres(npts)=R(npts)-data->domainLength;
  898. //Pres(npts)=Rdot(npts);
  899. }
  900. /*******************************************************************/
  901. //for (j = 1; j <=npts; j++) {
  902. // //for (k = 1; k <=nsp; k++) {
  903. // // Yres(j,k)=Ydot(j,k);
  904. // //}
  905. // //Tres(j)=Tdot(j);
  906. //}
  907. return(0);
  908. }
  909. void printSpaceTimeHeader(UserData data)
  910. {
  911. fprintf((data->output), "%15s\t","#1");
  912. for (size_t k = 1; k <=data->nvar+2; k++) {
  913. fprintf((data->output), "%15lu\t",k+1);
  914. }
  915. fprintf((data->output), "\n");
  916. fprintf((data->output), "%15s\t%15s\t%15s\t","#psi","time(s)","dpsi");
  917. fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
  918. for (size_t k = 1; k <=data->nsp; k++) {
  919. fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());
  920. }
  921. fprintf((data->output), "%15s\t","HRR(J)");
  922. fprintf((data->output), "%15s\n","Pressure(Pa)");
  923. }
  924. void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
  925. {
  926. double *ydata,*psidata;
  927. ydata = N_VGetArrayPointer_OpenMP(*y);
  928. if(data->adaptiveGrid){
  929. psidata = data->grid->x;
  930. }else{
  931. psidata = data->uniformGrid;
  932. }
  933. for (size_t i = 0; i < data->npts; i++) {
  934. fprintf(output, "%15.6e\t%15.6e\t",psi(i+1),t);
  935. if(i==0){
  936. fprintf(output, "%15.6e\t",psi(2)-psi(1));
  937. }
  938. else{
  939. fprintf(output, "%15.6e\t",psi(i+1)-psi(i));
  940. }
  941. for (size_t j = 0; j < data->nvar; j++) {
  942. fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
  943. }
  944. fprintf(output,"%15.6e\t",data->HRRdata[i]);
  945. fprintf(output, "\n");
  946. }
  947. fprintf(output, "\n\n");
  948. }
  949. void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data){
  950. double *ydata,*psidata, *ydotdata;
  951. ydata = N_VGetArrayPointer_OpenMP(*y);
  952. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  953. if(data->adaptiveGrid){
  954. psidata = data->grid->x;
  955. }else{
  956. psidata = data->uniformGrid;
  957. }
  958. fwrite(&t, sizeof(t), 1, output); //write time
  959. fwrite(psidata, data->npts*sizeof(psidata), 1, output); //write grid
  960. fwrite(ydata, data->neq*sizeof(ydata), 1, output); //write solution
  961. fwrite(ydotdata, data->neq*sizeof(ydotdata), 1, output); //write solutiondot
  962. }
  963. void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data){
  964. double *ydata,*psidata, *ydotdata;
  965. double t;
  966. if(data->adaptiveGrid){
  967. psidata = data->grid->x;
  968. }else{
  969. psidata = data->uniformGrid;
  970. }
  971. ydata = N_VGetArrayPointer_OpenMP(*y);
  972. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  973. fread(&t, sizeof(t), 1, input);
  974. data->tNow=t;
  975. fread(psidata, data->npts*sizeof(psidata), 1, input);
  976. fread(ydata, data->neq*sizeof(ydata), 1, input);
  977. fread(ydotdata, data->neq*sizeof(ydotdata), 1, input);
  978. if(data->adaptiveGrid){
  979. storeGrid(data->grid->x,data->grid->xOld,data->npts);
  980. }
  981. }
  982. void printGlobalHeader(UserData data)
  983. {
  984. fprintf((data->globalOutput), "%8s\t","#Time(s)");
  985. //fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");
  986. fprintf((data->globalOutput), "%15s","bFlux(kg/m^2/s)");
  987. fprintf((data->globalOutput), "%16s"," IsothermPos(m)");
  988. fprintf((data->globalOutput), "%15s\t","Pressure(Pa)");
  989. fprintf((data->globalOutput), "%15s\t","Pdot(Pa/s)");
  990. fprintf((data->globalOutput), "%15s\t","gamma");
  991. fprintf((data->globalOutput), "%15s\t","S_u(m/s)");
  992. fprintf((data->globalOutput), "%15s\t","Tu(K)");
  993. fprintf((data->globalOutput), "\n");
  994. }
  995. //
  996. //void printSpaceTimeRates(double t, N_Vector ydot, UserData data)
  997. //{
  998. // double *ydotdata,*psidata;
  999. // ydotdata = N_VGetArrayPointer_OpenMP(ydot);
  1000. // psidata = N_VGetArrayPointer_OpenMP(data->grid);
  1001. // for (int i = 0; i < data->npts; i++) {
  1002. // fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
  1003. // for (int j = 0; j < data->nvar; j++) {
  1004. // fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);
  1005. // }
  1006. // fprintf((data->ratesOutput), "\n");
  1007. // }
  1008. // fprintf((data->ratesOutput), "\n\n");
  1009. //}
  1010. //
  1011. void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)
  1012. {
  1013. double *ydata,*ydotdata, *innerMassFractionsData, *psidata;
  1014. innerMassFractionsData = data->innerMassFractions;
  1015. if(data->adaptiveGrid){
  1016. psidata = data->grid->x;
  1017. }else{
  1018. psidata = data->uniformGrid;
  1019. }
  1020. double TAvg, RAvg, YAvg, psiAvg;
  1021. ydata = N_VGetArrayPointer_OpenMP(*y);
  1022. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  1023. TAvg=data->isotherm;
  1024. double sum=ZERO;
  1025. double dpsim,area,aream,drdt;
  1026. double Cpb,Cvb,gamma,rho,flameArea,Tu;
  1027. /*Find the isotherm chosen by the user*/
  1028. size_t j=1;
  1029. size_t jj=1;
  1030. size_t jjj=1;
  1031. double wdot[data->nsp];
  1032. double wdotMax=0.0e0;
  1033. double advTerm=0.0e0;
  1034. psiAvg=0.0e0;
  1035. if(T(data->npts)>T(1)){
  1036. while (T(j)<TAvg) {
  1037. j=j+1;
  1038. }
  1039. YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(data->npts,data->k_oxidizer);
  1040. while (fabs((T(jj+1)-T(jj))/T(jj))>1e-08) {
  1041. jj=jj+1;
  1042. }
  1043. setGas(data,ydata,jj);
  1044. Tu=T(jj);
  1045. rho=data->gas->density();
  1046. Cpb=data->gas->cp_mass(); //J/kg/K
  1047. Cvb=data->gas->cv_mass(); //J/kg/K
  1048. gamma=Cpb/Cvb;
  1049. }
  1050. else{
  1051. while (T(j)>TAvg) {
  1052. j=j+1;
  1053. }
  1054. YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(1,data->k_oxidizer);
  1055. while (fabs((T(data->npts-jj-1)-T(data->npts-jj))/T(data->npts-jj))>1e-08) {
  1056. jj=jj+1;
  1057. }
  1058. setGas(data,ydata,data->npts-jj);
  1059. Tu=T(data->npts-jj);
  1060. rho=data->gas->density();
  1061. Cpb=data->gas->cp_mass(); //J/kg/K
  1062. Cvb=data->gas->cv_mass(); //J/kg/K
  1063. gamma=Cpb/Cvb;
  1064. }
  1065. if(T(j)<TAvg){
  1066. RAvg=((R(j+1)-R(j))/(T(j+1)-T(j)))*(TAvg-T(j))+R(j);
  1067. }
  1068. else{
  1069. RAvg=((R(j)-R(j-1))/(T(j)-T(j-1)))*(TAvg-T(j-1))+R(j-1);
  1070. }
  1071. ////Experimental burning speed calculation:
  1072. //int nMax=0;
  1073. ////nMax=maxCurvIndex(ydata, data->nt, data->nvar,
  1074. //// data->grid->x, data->npts);
  1075. //nMax=maxGradIndex(ydata, data->nt, data->nvar,
  1076. // data->grid->x, data->npts);
  1077. //advTerm=(T(nMax)-T(nMax-1))/(data->mass*(psi(nMax)-psi(nMax-1)));
  1078. //aream=calc_area(R(nMax),&data->metric);
  1079. ////setGas(data,ydata,nMax);
  1080. ////rho=data->gas->density();
  1081. //psiAvg=-Tdot(nMax)/(rho*aream*advTerm);
  1082. ////if(t>data->ignTime){
  1083. //// for(size_t n=2;n<data->npts;n++){
  1084. //// setGas(data,ydata,n);
  1085. //// data->gas->getNetProductionRates(wdot); //kmol/m^3
  1086. //// advTerm=(T(n)-T(n-1))/(data->mass*(psi(n)-psi(n-1)));
  1087. //// if(fabs(wdot[data->k_oxidizer-1])>=wdotMax){
  1088. //// aream=calc_area(R(n),&data->metric);
  1089. //// psiAvg=-Tdot(n)/(rho*aream*advTerm);
  1090. //// wdotMax=fabs(wdot[data->k_oxidizer-1]);
  1091. //// }
  1092. //// }
  1093. ////}
  1094. ////else{
  1095. //// psiAvg=0.0e0;
  1096. ////}
  1097. //drdt=(RAvg-data->flamePosition[1])/(t-data->flameTime[1]);
  1098. //data->flamePosition[0]=data->flamePosition[1];
  1099. //data->flamePosition[1]=RAvg;
  1100. //data->flameTime[0]=data->flameTime[1];
  1101. //data->flameTime[1]=t;
  1102. //flameArea=calc_area(RAvg,&data->metric);
  1103. /*Use the Trapezoidal rule to calculate the mass burning rate based on
  1104. * the consumption of O2*/
  1105. aream= calc_area(R(1)+1e-03*data->domainLength,&data->metric);
  1106. for (j = 2; j <data->npts; j++) {
  1107. dpsim=(psi(j)-psi(j-1))*data->mass;
  1108. area= calc_area(R(j),&data->metric);
  1109. sum=sum+HALF*dpsim*((Ydot(j-1,data->k_oxidizer)/aream)
  1110. +(Ydot(j,data->k_oxidizer)/area));
  1111. aream=area;
  1112. }
  1113. //double maxOH,maxHO2;
  1114. //maxOH=0.0e0;
  1115. //maxHO2=0.0e0;
  1116. //for(j=1;j<data->npts;j++){
  1117. // if(Y(j,data->k_OH)>maxOH){
  1118. // maxOH=Y(j,data->k_OH);
  1119. // }
  1120. //}
  1121. //for(j=1;j<data->npts;j++){
  1122. // if(Y(j,data->k_HO2)>maxHO2){
  1123. // maxHO2=Y(j,data->k_HO2);
  1124. // }
  1125. //}
  1126. fprintf((data->globalOutput), "%15.6e\t",t);
  1127. //fprintf((data->globalOutput), "%15.6e\t",psiAvg);
  1128. fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/YAvg);
  1129. fprintf((data->globalOutput), "%15.6e\t",RAvg);
  1130. fprintf((data->globalOutput), "%15.6e\t",P(data->npts));
  1131. fprintf((data->globalOutput), "%15.6e\t",Pdot(data->npts));
  1132. fprintf((data->globalOutput), "%15.6e\t",gamma);
  1133. fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/(YAvg*rho));
  1134. fprintf((data->globalOutput), "%15.6e\t",Tu);
  1135. fprintf((data->globalOutput), "\n");
  1136. }
  1137. //
  1138. //void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data)
  1139. //{
  1140. // double *ydata,*psidata;
  1141. // ydata = N_VGetArrayPointer_OpenMP(y);
  1142. // psidata = N_VGetArrayPointer_OpenMP(data->grid);
  1143. // for (int i = 0; i < data->npts; i++) {
  1144. // fprintf((data->gridOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
  1145. // for (int j = 0; j < data->nvar; j++) {
  1146. // fprintf((data->gridOutput), "%15.6e\t",ydata[j+i*data->nvar]);
  1147. // }
  1148. // fprintf((data->gridOutput), "\n");
  1149. // }
  1150. // fprintf((data->gridOutput), "\n\n");
  1151. //}
  1152. //
  1153. //
  1154. ////void repairSolution(N_Vector y, N_Vector ydot, UserData data){
  1155. //// int npts=data->npts;
  1156. //// double *ydata;
  1157. //// double *ydotdata;
  1158. //// ydata = N_VGetArrayPointer_OpenMP(y);
  1159. //// ydotdata = N_VGetArrayPointer_OpenMP(ydot);
  1160. ////
  1161. //// T(2)=T(1);
  1162. //// T(npts-1)=T(npts);
  1163. //// Tdot(2)=Tdot(1);
  1164. //// Tdot(npts-1)=Tdot(npts);
  1165. //// for (int k = 1; k <=data->nsp; k++) {
  1166. //// Y(2,k)=Y(1,k);
  1167. //// Y(npts-1,k)=Y(npts,k);
  1168. ////
  1169. //// Ydot(2,k)=Ydot(1,k);
  1170. //// Ydot(npts-1,k)=Ydot(npts,k);
  1171. //// }
  1172. ////}