A low Mach, 1D, reacting flow code.
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

1303 linhas
37KB

  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. dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;
  660. double mass, mdotIn;
  661. double sum, sum1, sum2, sum3;
  662. size_t j,k;
  663. int m;
  664. m=data->metric; //Unitless
  665. mass=data->mass; //Units: kg
  666. mdotIn=data->mdot*calc_area(R(npts),&m); //Units: kg/s
  667. // /*evaluate properties at j=1*************************/
  668. setGas(data,ydata,1);
  669. rhom=data->gas->density();
  670. Cpb=data->gas->cp_mass(); //J/kg/K
  671. Cvb=data->gas->cv_mass(); //J/kg/K
  672. aream= calc_area(R(1),&m);
  673. /*******************************************************************/
  674. /*Calculate values at j=2's m and mhalf*****************************/
  675. getTransport(data, ydata, 1, &rhomhalf,&lambdamhalf,YVmhalf);
  676. areamhalf= calc_area(HALF*(R(1)+R(2)),&m);
  677. areamhalfsq= areamhalf*areamhalf;
  678. /*******************************************************************/
  679. /*Fill up res with left side (center) boundary conditions:**********/
  680. /*We impose zero fluxes at the center:*/
  681. /*Mass:*/
  682. Rres(1)=Rdot(1);
  683. /*Energy:*/
  684. if(data->dirichletInner){
  685. //Tres(1)=Tdot(1);
  686. Tres(1)=T(1)-data->innerTemperature;
  687. }
  688. else{
  689. Tres(1)=T(2)-T(1);
  690. //Tres(1)=Tdot(1) - (Pdot(1)/(rhom*Cpb))
  691. // +(double)(data->metric+1)*(rhomhalf*lambdamhalf*areamhalfsq*(T(2)-T(1))/psi(2)-psi(1));
  692. }
  693. /*Species:*/
  694. sum=ZERO;
  695. for (k = 1; k <=nsp; k++) {
  696. if(k!=k_bath){
  697. if(fabs(mdotIn)>1e-14){
  698. Yres(1,k)=innerMassFractionsData[k-1]-
  699. Y(1,k)-
  700. (YVmhalf(k)*areamhalf)/mdotIn;
  701. }
  702. else{
  703. //Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];
  704. Yres(1,k)=Y(2,k)-Y(1,k);
  705. /*The below flux boundary condition makes the
  706. * problem more prone to diverge. How does one
  707. * fix this?*/
  708. //Yres(1,k)=YVmhalf(k);
  709. }
  710. sum=sum+Y(1,k);
  711. }
  712. }
  713. Yres(1,k_bath)=ONE-sum-Y(1,k_bath);
  714. /*Pressure:*/
  715. Pres(1)=P(2)-P(1);
  716. /*Fill up res with governing equations at inner points:*************/
  717. for (j = 2; j < npts; j++) {
  718. /*evaluate various mesh differences*///
  719. dpsip = (psi(j+1) - psi(j) )*mass;
  720. dpsim = (psi(j) - psi(j-1))*mass;
  721. dpsiav = HALF*(psi(j+1) - psi(j-1))*mass;
  722. dpsipm = (psi(j+1) - psi(j-1))*mass;
  723. /***********************************///
  724. /*evaluate various central difference coefficients*/
  725. cendfm = - dpsip / (dpsim*dpsipm);
  726. cendfc = (dpsip-dpsim) / (dpsip*dpsim);
  727. cendfp = dpsim / (dpsip*dpsipm);
  728. /**************************************************/
  729. /*evaluate properties at j*************************/
  730. setGas(data,ydata,j);
  731. rho=data->gas->density(); //kg/m^3
  732. Cpb=data->gas->cp_mass(); //J/kg/K
  733. Cvb=data->gas->cv_mass(); //J/kg/K
  734. data->gas->getNetProductionRates(wdot); //kmol/m^3
  735. data->gas->getEnthalpy_RT(enthalpy); //unitless
  736. data->gas->getCp_R(Cp); //unitless
  737. area = calc_area(R(j),&m); //m^2
  738. /*evaluate properties at p*************************/
  739. getTransport(data, ydata, j, &rhophalf,&lambdaphalf,YVphalf);
  740. areaphalf= calc_area(HALF*(R(j)+R(j+1)),&m);
  741. areaphalfsq= areaphalf*areaphalf;
  742. /**************************************************///
  743. /*Mass:*/
  744. /* ∂r/∂ψ = 1/ρA */
  745. Rres(j)=((R(j)-R(j-1))/dpsim)-(TWO/(rhom*aream+rho*area));
  746. /*Energy:*/
  747. /* ∂T/∂t = - ṁ(∂T/∂ψ)
  748. * + (1/cₚ)(∂/∂ψ)(λρA²∂T/∂ψ)
  749. * - (A/cₚ) ∑ YᵢVᵢcₚᵢ(∂T/∂ψ)
  750. * - (1/ρcₚ)∑ ώᵢhᵢ
  751. * + (1/ρcₚ)(∂P/∂t) */
  752. /*Notes:
  753. * λ has units J/m/s/K.
  754. * YᵢVᵢ has units kg/m^2/s.
  755. * hᵢ has units J/kmol, so we must multiply the enthalpy
  756. * defined above (getEnthalpy_RT) by T (K) and the gas constant
  757. * (J/kmol/K) to get the right units.
  758. * cₚᵢ has units J/kg/K, so we must multiply the specific heat
  759. * defined above (getCp_R) by the gas constant (J/kmol/K) and
  760. * divide by the molecular weight (kg/kmol) to get the right
  761. * units.
  762. * */
  763. //enthalpy formulation:
  764. tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));
  765. sum=ZERO;
  766. sum1=ZERO;
  767. for (k = 1; k <=nsp; k++) {
  768. sum=sum+wdot(k)*enthalpy(k);
  769. sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*HALF*(YVmhalf(k)+YVphalf(k));
  770. }
  771. sum=sum*Cantera::GasConstant*T(j);
  772. sum1=sum1*Cantera::GasConstant;
  773. diffTerm =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
  774. -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
  775. /(dpsiav*Cpb) )
  776. -(sum1*area*(cendfp*T(j+1)
  777. +cendfc*T(j)
  778. +cendfm*T(j-1))/Cpb);
  779. srcTerm = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cpb);
  780. advTerm = (mdotIn*(T(j)-T(j-1))/dpsim);
  781. Tres(j)= tranTerm
  782. +advTerm
  783. -diffTerm
  784. +srcTerm;
  785. // //energy formulation:
  786. // tranTerm = Tdot(j);
  787. // sum=ZERO;
  788. // sum1=ZERO;
  789. // sum2=ZERO;
  790. // sum3=ZERO;
  791. // for (k = 1; k <=nsp; k++) {
  792. // energy(k)=enthalpy(k)-ONE;
  793. // sum=sum+wdot(k)*energy(k);
  794. // sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
  795. // *HALF*(YVmhalf(k)+YVphalf(k));
  796. // sum2=sum2+(YVmhalf(k)/data->gas->molecularWeight(k-1));
  797. // sum3=sum3+(YVphalf(k)/data->gas->molecularWeight(k-1));
  798. // }
  799. // sum=sum*Cantera::GasConstant*T(j);
  800. // sum1=sum1*Cantera::GasConstant;
  801. // diffTerm =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
  802. // -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
  803. // /(dpsiav*Cvb) )
  804. // -(sum1*area*(cendfp*T(j+1)
  805. // +cendfc*T(j)
  806. // +cendfm*T(j-1))/Cvb);
  807. // srcTerm = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cvb);
  808. // advTerm = (mdotIn*(T(j)-T(j-1))/dpsim);
  809. // advTerm = advTerm + (Cantera::GasConstant*T(j)*area/Cvb)*((sum3-sum2)/dpsiav);
  810. // Tres(j)= tranTerm
  811. // +advTerm
  812. // -diffTerm
  813. // +srcTerm;
  814. /*Species:*/
  815. /* ∂Yᵢ/∂t = - ṁ(∂Yᵢ/∂ψ)
  816. * - (∂/∂ψ)(AYᵢVᵢ)
  817. * + (ώᵢWᵢ/ρ) */
  818. sum=ZERO;
  819. for (k = 1; k <=nsp; k++) {
  820. if(k!=k_bath){
  821. tranTerm = Ydot(j,k);
  822. diffTerm = (YVphalf(k)*areaphalf
  823. -YVmhalf(k)*areamhalf)/dpsiav;
  824. srcTerm = wdot(k)
  825. *(data->gas->molecularWeight(k-1))/rho;
  826. advTerm = (mdotIn*(Y(j,k)-Y(j-1,k))/dpsim);
  827. Yres(j,k)= tranTerm
  828. +advTerm
  829. +diffTerm
  830. -srcTerm;
  831. sum=sum+Y(j,k);
  832. }
  833. }
  834. Yres(j,k_bath)=ONE-sum-Y(j,k_bath);
  835. /*Pressure:*/
  836. Pres(j) = P(j+1)-P(j);
  837. /*Assign values evaluated at p and phalf to m
  838. * and mhalf to save some cpu cost:****************/
  839. areamhalf=areaphalf;
  840. areamhalfsq=areaphalfsq;
  841. aream=area;
  842. rhom=rho;
  843. rhomhalf=rhophalf;
  844. lambdamhalf=lambdaphalf;
  845. for (k = 1; k <=nsp; k++) {
  846. YVmhalf(k)=YVphalf(k);
  847. }
  848. /**************************************************/
  849. }
  850. /*******************************************************************///
  851. /*Fill up res with right side (wall) boundary conditions:***********/
  852. /*We impose zero fluxes at the wall:*/
  853. setGas(data,ydata,npts);
  854. rho=data->gas->density();
  855. area = calc_area(R(npts),&m);
  856. /*Mass:*/
  857. dpsim=(psi(npts)-psi(npts-1))*mass;
  858. Rres(npts)=((R(npts)-R(npts-1))/dpsim)-(TWO/(rhom*aream+rho*area));
  859. /*Energy:*/
  860. if(data->dirichletOuter){
  861. Tres(npts)=T(npts)-data->wallTemperature;
  862. }
  863. else{
  864. Tres(npts)=T(npts)-T(npts-1);
  865. }
  866. /*Species:*/
  867. sum=ZERO;
  868. for (k = 1; k <=nsp; k++) {
  869. if(k!=k_bath){
  870. //Yres(npts,k)=Y(npts,k)-Y(npts-1,k);
  871. Yres(npts,k)=YVmhalf(k);
  872. sum=sum+Y(npts,k);
  873. }
  874. }
  875. Yres(npts,k_bath)=ONE-sum-Y(npts,k_bath);
  876. /*Pressure:*/
  877. if(data->constantPressure){
  878. Pres(npts)=Pdot(npts)-data->dPdt;
  879. }
  880. else{
  881. Pres(npts)=R(npts)-data->domainLength;
  882. //Pres(npts)=Rdot(npts);
  883. }
  884. /*******************************************************************/
  885. //for (j = 1; j <=npts; j++) {
  886. // //for (k = 1; k <=nsp; k++) {
  887. // // Yres(j,k)=Ydot(j,k);
  888. // //}
  889. // //Tres(j)=Tdot(j);
  890. //}
  891. return(0);
  892. }
  893. void printSpaceTimeHeader(UserData data)
  894. {
  895. fprintf((data->output), "%15s\t","#1");
  896. for (size_t k = 1; k <=data->nvar+2; k++) {
  897. fprintf((data->output), "%15lu\t",k+1);
  898. }
  899. fprintf((data->output), "\n");
  900. fprintf((data->output), "%15s\t%15s\t%15s\t","#psi","time(s)","dpsi");
  901. fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
  902. for (size_t k = 1; k <=data->nsp; k++) {
  903. fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());
  904. }
  905. fprintf((data->output), "%15s\n","Pressure(Pa)");
  906. }
  907. void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
  908. {
  909. double *ydata,*psidata;
  910. ydata = N_VGetArrayPointer_OpenMP(*y);
  911. if(data->adaptiveGrid){
  912. psidata = data->grid->x;
  913. }else{
  914. psidata = data->uniformGrid;
  915. }
  916. for (size_t i = 0; i < data->npts; i++) {
  917. fprintf(output, "%15.6e\t%15.6e\t",psi(i+1),t);
  918. if(i==0){
  919. fprintf(output, "%15.6e\t",psi(2)-psi(1));
  920. }
  921. else{
  922. fprintf(output, "%15.6e\t",psi(i+1)-psi(i));
  923. }
  924. for (size_t j = 0; j < data->nvar; j++) {
  925. fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
  926. }
  927. fprintf(output, "\n");
  928. }
  929. fprintf(output, "\n\n");
  930. }
  931. void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data){
  932. double *ydata,*psidata, *ydotdata;
  933. ydata = N_VGetArrayPointer_OpenMP(*y);
  934. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  935. if(data->adaptiveGrid){
  936. psidata = data->grid->x;
  937. }else{
  938. psidata = data->uniformGrid;
  939. }
  940. fwrite(&t, sizeof(t), 1, output); //write time
  941. fwrite(psidata, data->npts*sizeof(psidata), 1, output); //write grid
  942. fwrite(ydata, data->neq*sizeof(ydata), 1, output); //write solution
  943. fwrite(ydotdata, data->neq*sizeof(ydotdata), 1, output); //write solutiondot
  944. }
  945. void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data){
  946. double *ydata,*psidata, *ydotdata;
  947. double t;
  948. if(data->adaptiveGrid){
  949. psidata = data->grid->x;
  950. }else{
  951. psidata = data->uniformGrid;
  952. }
  953. ydata = N_VGetArrayPointer_OpenMP(*y);
  954. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  955. fread(&t, sizeof(t), 1, input);
  956. data->tNow=t;
  957. fread(psidata, data->npts*sizeof(psidata), 1, input);
  958. fread(ydata, data->neq*sizeof(ydata), 1, input);
  959. fread(ydotdata, data->neq*sizeof(ydotdata), 1, input);
  960. if(data->adaptiveGrid){
  961. storeGrid(data->grid->x,data->grid->xOld,data->npts);
  962. }
  963. }
  964. void printGlobalHeader(UserData data)
  965. {
  966. fprintf((data->globalOutput), "%8s\t","#Time(s)");
  967. //fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");
  968. fprintf((data->globalOutput), "%15s","bFlux(kg/m^2/s)");
  969. fprintf((data->globalOutput), "%16s"," IsothermPos(m)");
  970. fprintf((data->globalOutput), "%15s\t","Pressure(Pa)");
  971. fprintf((data->globalOutput), "%15s\t","Pdot(Pa/s)");
  972. fprintf((data->globalOutput), "%15s\t","gamma");
  973. fprintf((data->globalOutput), "%15s\t","S_u(m/s)");
  974. fprintf((data->globalOutput), "%15s\t","Tu(K)");
  975. fprintf((data->globalOutput), "\n");
  976. }
  977. //
  978. //void printSpaceTimeRates(double t, N_Vector ydot, UserData data)
  979. //{
  980. // double *ydotdata,*psidata;
  981. // ydotdata = N_VGetArrayPointer_OpenMP(ydot);
  982. // psidata = N_VGetArrayPointer_OpenMP(data->grid);
  983. // for (int i = 0; i < data->npts; i++) {
  984. // fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
  985. // for (int j = 0; j < data->nvar; j++) {
  986. // fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);
  987. // }
  988. // fprintf((data->ratesOutput), "\n");
  989. // }
  990. // fprintf((data->ratesOutput), "\n\n");
  991. //}
  992. //
  993. void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)
  994. {
  995. double *ydata,*ydotdata, *innerMassFractionsData, *psidata;
  996. innerMassFractionsData = data->innerMassFractions;
  997. if(data->adaptiveGrid){
  998. psidata = data->grid->x;
  999. }else{
  1000. psidata = data->uniformGrid;
  1001. }
  1002. double TAvg, RAvg, YAvg, psiAvg;
  1003. ydata = N_VGetArrayPointer_OpenMP(*y);
  1004. ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
  1005. TAvg=data->isotherm;
  1006. double sum=ZERO;
  1007. double dpsim,area,aream,drdt;
  1008. double Cpb,Cvb,gamma,rho,flameArea,Tu;
  1009. /*Find the isotherm chosen by the user*/
  1010. size_t j=1;
  1011. size_t jj=1;
  1012. size_t jjj=1;
  1013. double wdot[data->nsp];
  1014. double wdotMax=0.0e0;
  1015. double advTerm=0.0e0;
  1016. psiAvg=0.0e0;
  1017. if(T(data->npts)>T(1)){
  1018. while (T(j)<TAvg) {
  1019. j=j+1;
  1020. }
  1021. YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(data->npts,data->k_oxidizer);
  1022. while (fabs((T(jj+1)-T(jj))/T(jj))>1e-08) {
  1023. jj=jj+1;
  1024. }
  1025. setGas(data,ydata,jj);
  1026. Tu=T(jj);
  1027. rho=data->gas->density();
  1028. Cpb=data->gas->cp_mass(); //J/kg/K
  1029. Cvb=data->gas->cv_mass(); //J/kg/K
  1030. gamma=Cpb/Cvb;
  1031. }
  1032. else{
  1033. while (T(j)>TAvg) {
  1034. j=j+1;
  1035. }
  1036. YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(1,data->k_oxidizer);
  1037. while (fabs((T(data->npts-jj-1)-T(data->npts-jj))/T(data->npts-jj))>1e-08) {
  1038. jj=jj+1;
  1039. }
  1040. setGas(data,ydata,data->npts-jj);
  1041. Tu=T(data->npts-jj);
  1042. rho=data->gas->density();
  1043. Cpb=data->gas->cp_mass(); //J/kg/K
  1044. Cvb=data->gas->cv_mass(); //J/kg/K
  1045. gamma=Cpb/Cvb;
  1046. }
  1047. if(T(j)<TAvg){
  1048. RAvg=((R(j+1)-R(j))/(T(j+1)-T(j)))*(TAvg-T(j))+R(j);
  1049. }
  1050. else{
  1051. RAvg=((R(j)-R(j-1))/(T(j)-T(j-1)))*(TAvg-T(j-1))+R(j-1);
  1052. }
  1053. ////Experimental burning speed calculation:
  1054. //int nMax=0;
  1055. ////nMax=maxCurvIndex(ydata, data->nt, data->nvar,
  1056. //// data->grid->x, data->npts);
  1057. //nMax=maxGradIndex(ydata, data->nt, data->nvar,
  1058. // data->grid->x, data->npts);
  1059. //advTerm=(T(nMax)-T(nMax-1))/(data->mass*(psi(nMax)-psi(nMax-1)));
  1060. //aream=calc_area(R(nMax),&data->metric);
  1061. ////setGas(data,ydata,nMax);
  1062. ////rho=data->gas->density();
  1063. //psiAvg=-Tdot(nMax)/(rho*aream*advTerm);
  1064. ////if(t>data->ignTime){
  1065. //// for(size_t n=2;n<data->npts;n++){
  1066. //// setGas(data,ydata,n);
  1067. //// data->gas->getNetProductionRates(wdot); //kmol/m^3
  1068. //// advTerm=(T(n)-T(n-1))/(data->mass*(psi(n)-psi(n-1)));
  1069. //// if(fabs(wdot[data->k_oxidizer-1])>=wdotMax){
  1070. //// aream=calc_area(R(n),&data->metric);
  1071. //// psiAvg=-Tdot(n)/(rho*aream*advTerm);
  1072. //// wdotMax=fabs(wdot[data->k_oxidizer-1]);
  1073. //// }
  1074. //// }
  1075. ////}
  1076. ////else{
  1077. //// psiAvg=0.0e0;
  1078. ////}
  1079. //drdt=(RAvg-data->flamePosition[1])/(t-data->flameTime[1]);
  1080. //data->flamePosition[0]=data->flamePosition[1];
  1081. //data->flamePosition[1]=RAvg;
  1082. //data->flameTime[0]=data->flameTime[1];
  1083. //data->flameTime[1]=t;
  1084. //flameArea=calc_area(RAvg,&data->metric);
  1085. /*Use the Trapezoidal rule to calculate the mass burning rate based on
  1086. * the consumption of O2*/
  1087. aream= calc_area(R(1)+1e-03*data->domainLength,&data->metric);
  1088. for (j = 2; j <data->npts; j++) {
  1089. dpsim=(psi(j)-psi(j-1))*data->mass;
  1090. area= calc_area(R(j),&data->metric);
  1091. sum=sum+HALF*dpsim*((Ydot(j-1,data->k_oxidizer)/aream)
  1092. +(Ydot(j,data->k_oxidizer)/area));
  1093. aream=area;
  1094. }
  1095. //double maxOH,maxHO2;
  1096. //maxOH=0.0e0;
  1097. //maxHO2=0.0e0;
  1098. //for(j=1;j<data->npts;j++){
  1099. // if(Y(j,data->k_OH)>maxOH){
  1100. // maxOH=Y(j,data->k_OH);
  1101. // }
  1102. //}
  1103. //for(j=1;j<data->npts;j++){
  1104. // if(Y(j,data->k_HO2)>maxHO2){
  1105. // maxHO2=Y(j,data->k_HO2);
  1106. // }
  1107. //}
  1108. fprintf((data->globalOutput), "%15.6e\t",t);
  1109. //fprintf((data->globalOutput), "%15.6e\t",psiAvg);
  1110. fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/YAvg);
  1111. fprintf((data->globalOutput), "%15.6e\t",RAvg);
  1112. fprintf((data->globalOutput), "%15.6e\t",P(data->npts));
  1113. fprintf((data->globalOutput), "%15.6e\t",Pdot(data->npts));
  1114. fprintf((data->globalOutput), "%15.6e\t",gamma);
  1115. fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/(YAvg*rho));
  1116. fprintf((data->globalOutput), "%15.6e\t",Tu);
  1117. fprintf((data->globalOutput), "\n");
  1118. }
  1119. //
  1120. //void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data)
  1121. //{
  1122. // double *ydata,*psidata;
  1123. // ydata = N_VGetArrayPointer_OpenMP(y);
  1124. // psidata = N_VGetArrayPointer_OpenMP(data->grid);
  1125. // for (int i = 0; i < data->npts; i++) {
  1126. // fprintf((data->gridOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
  1127. // for (int j = 0; j < data->nvar; j++) {
  1128. // fprintf((data->gridOutput), "%15.6e\t",ydata[j+i*data->nvar]);
  1129. // }
  1130. // fprintf((data->gridOutput), "\n");
  1131. // }
  1132. // fprintf((data->gridOutput), "\n\n");
  1133. //}
  1134. //
  1135. //
  1136. ////void repairSolution(N_Vector y, N_Vector ydot, UserData data){
  1137. //// int npts=data->npts;
  1138. //// double *ydata;
  1139. //// double *ydotdata;
  1140. //// ydata = N_VGetArrayPointer_OpenMP(y);
  1141. //// ydotdata = N_VGetArrayPointer_OpenMP(ydot);
  1142. ////
  1143. //// T(2)=T(1);
  1144. //// T(npts-1)=T(npts);
  1145. //// Tdot(2)=Tdot(1);
  1146. //// Tdot(npts-1)=Tdot(npts);
  1147. //// for (int k = 1; k <=data->nsp; k++) {
  1148. //// Y(2,k)=Y(1,k);
  1149. //// Y(npts-1,k)=Y(npts,k);
  1150. ////
  1151. //// Ydot(2,k)=Ydot(1,k);
  1152. //// Ydot(npts-1,k)=Ydot(npts,k);
  1153. //// }
  1154. ////}