weiye
/
Droplet_LTORC
关注 1
点赞 0
派生 0
代码 工单 0 合并请求 0 版本发布 0 百科 动态
Droplet Lagrangian Transient One-dimensional Reacting Code Implementation of both liquid and gas phase governing equations.
您最多选择25个主题 主题必须以字母或数字开头，可以包含连字符 (-)，并且长度不得超过35个字符
12 提交
1 分支
1.6MB
分支: binary_fuel
分支列表 标签列表
${ item.name }
创建分支 ${ searchTerm }
从 'binary_fuel'
${ noResults }
Droplet_LTORC/src/UserData.cpp

674 行
20KB
原始文件 永久链接 Blame 文件历史 

  1. #include "UserData.h"

  2. #include "parse.h"

  3. 

  4. void freeUserData(UserData data){

  5. 

  6. 	if(data!=NULL){

  7.         /*destructor for C++ vector*/

  8.         using Td = std::vector<double>;

  9.         using Ts = std::vector<std::string>;

  10.         using Tt = std::vector<size_t>;

  11.         data->dropMole.~Td();

  12.         data->dropSpec.~Ts();

  13.         data->gamma.~Td();

  14.         data->MW.~Td();

  15.         data->k_drop.~Tt();

  16. 

  17. 		if(data->trmix!=NULL){

  18. 			delete data->trmix;

  19. 			printf("Transport Deleted!\n");

  20. 

  21. 		}

  22. 		if(data->gas!=NULL){

  23. 			delete data->gas;

  24. 			printf("Gas Deleted!\n");

  25. 

  26. 		}

  27. 		if(data->adaptiveGrid){

  28. 			if(data->grid->xOld!=NULL){

  29. 				delete[] data->grid->xOld;

  30. 				printf("old grid array Deleted!\n");

  31. 			}

  32. 			if(data->grid->x!=NULL){

  33. 				delete[] data->grid->x;

  34. 				printf("current grid array Deleted!\n");

  35. 			}

  36. 			if(data->grid!=NULL){

  37. 				free(data->grid);

  38. 				printf("grid object Freed!\n");

  39. 			}

  40. 		}

  41. 		else{

  42. 			if(data->uniformGrid!=NULL){

  43. 				delete[] data->uniformGrid;

  44. 				printf("uniformGrid deleted!\n");

  45. 			}

  46. 		}

  47. 		if(data->innerMassFractions!=NULL){

  48. 			delete[] data->innerMassFractions;

  49. 			printf("innerMassFractions array Deleted!\n");

  50. 		}

  51. 		if(data->output!=NULL){

  52.   			fclose(data->output);

  53. 			printf("Output File Cleared from Memory!\n");

  54. 		}

  55. 		if(data->gridOutput!=NULL){

  56.   			fclose(data->gridOutput);

  57. 			printf("Grid Output File Cleared from Memory!\n");

  58. 		}

  59. 		//if(data->ratesOutput!=NULL){

  60.   		//	fclose(data->ratesOutput);

  61. 		//	printf("Rates Output File Cleared from Memory!\n");

  62. 		//}

  63. 		if(data->globalOutput!=NULL){

  64.   			fclose(data->globalOutput);

  65. 			printf("Global Output File Cleared from Memory!\n");

  66. 		}

  67. 		// if(data->timescaleOutput!=NULL){

  68.   		// 	fclose(data->timescaleOutput);

  69. 		// 	printf("Characteristic Timescale Output File Cleared from Memory!\n");

  70. 		// }

  71. 		//if(data->rxnROPOutput!=NULL){

  72.   		//	fclose(data->rxnROPOutput);

  73. 		//	printf("Reactions Rate of Progress Output File Cleared from Memory!\n");

  74. 		//}

  75. 		//if(data->spROPOutput!=NULL){

  76.   		//	fclose(data->spROPOutput);

  77. 		//	printf("Species Rate of Production Output File Cleared from Memory!\n");

  78. 		//}

  79. 	}

  80.   	free(data);             /* Free the user data */

  81. 	printf("\n\n");

  82. }

  83. 

  84. UserData allocateUserData(FILE *input){

  85. 

  86. 	UserData data;

  87.   	data = (UserData) malloc(sizeof *data);

  88.     new(&(data->dropMole)) std::vector<double>;

  89.     new(&(data->dropSpec)) std::vector<std::string>;

  90.     new(&(data->gamma))   std::vector<double>;

  91.     new(&(data->MW))   std::vector<double>;

  92.     new(&(data->k_drop))   std::vector<size_t>;

  93. 

  94. 	if(!data){

  95. 		printf("Allocation Failed!\n");

  96. 		return(NULL);

  97. 	}

  98. 	setSaneDefaults(data);

  99. 

  100. 	int ier;

  101. 

  102. 	ier=parseNumber<size_t>(input, "liquidBasePts"	, MAXBUFLEN, &data->l_npts);

  103. 	if(ier==-1 || data->l_npts<=0){

  104. 		printf("Enter non-zero liquidbasePts!\n");

  105. 		return(NULL);

  106. 	}

  107. 	ier=parseNumber<size_t>(input, "gasBasePts"	, MAXBUFLEN, &data->g_npts);

  108. 	if(ier==-1 || data->g_npts<=0){

  109. 		printf("Enter non-zero gasBasePts!\n");

  110. 		return(NULL);

  111. 	}

  112. 

  113. 	ier=parseNumber<double>(input, "domainLength", MAXBUFLEN, &data->domainLength);

  114. 	if(ier==-1 || data->domainLength<=0.0e0){

  115. 		printf("domainLength error!\n");

  116. 		return(NULL);

  117. 	}

  118. 

  119. 	ier=parseNumber<double>(input, "Rd", MAXBUFLEN, &data->Rd);

  120. 	if(ier==-1 || data->Rd<=0.0e0){

  121. 		printf("Rd error!\n");

  122. 		return(NULL);

  123. 	}

  124. 

  125. 	ier=parseNumber<int>(input, "constantPressure"	, MAXBUFLEN, &data->constantPressure);

  126. 	if(ier==-1 || (data->constantPressure!=0 && data->constantPressure!=1)){

  127. 		printf("constantPressure error!\n");

  128. 		return(NULL);

  129. 	}

  130. 

  131. 	ier=parseNumber<int>(input, "problemType"	, MAXBUFLEN, &data->problemType);

  132. 	if(ier==-1 || (data->problemType!=0 

  133. 		       && data->problemType!=1

  134. 			&& data->problemType!=2

  135. 			&& data->problemType!=3)){

  136. 		printf("include valid problemType!\n");

  137. 		printf("0: premixed combustion with NO equilibrated ignition kernel\n");

  138. 		printf("1: premixed combustion WITH equilibrated ignition kernel\n");

  139. 		printf("2: arbitrary initial conditions\n");

  140. 		printf("3: Restart\n");

  141. 		return(NULL);

  142. 	}

  143. 

  144. 	ier=parseNumber<int>(input,"dropletType",MAXBUFLEN,&data->dropType);

  145. 	if(ier==-1 || (data->dropType!=0 

  146. 		       && data->dropType!=1)){

  147. 		printf("include valid dropletType!\n");

  148. 		printf("0: single component droplet\n");

  149. 		printf("1: bi-component droplet\n");

  150. 		return(NULL);

  151. 	}

  152. 

  153. 	ier=parseDropSpec<double>(input,"dropletMole",MAXBUFLEN,&data->dropType,data->dropMole);

  154. 	if(ier== -1){

  155. 		printf("include correct format for droplet mole fraction!\n");

  156. 		printf("single composition droplet: dropletMole:1.00\n");

  157. 		printf("binary composition droplet: dropletMole:0.3,0.7\n");

  158. 		return(NULL);

  159. 	}

  160. 

  161. 	ier=parseDropSpec<std::string>(input,"dropletComposition",MAXBUFLEN,&data->dropType,data->dropSpec);

  162. 	if(ier== -1){

  163. 		printf("include correct format for droplet composition!\n");

  164. 		printf("single composition droplet: dropletComposition:n-heptane\n");

  165. 		printf("binary composition droplet: dropletComposition:propane,n-heptane\n");

  166. 		return(NULL);

  167. 	}

  168. 	

  169. 	

  170. 	ier=parseNumber<int>(input, "quasiSteady"	, MAXBUFLEN, &data->quasiSteady);

  171. 	if(ier==-1 || (data->quasiSteady!=0 

  172. 		       && data->quasiSteady!=1)){

  173. 		printf("include valid quasiSteady!\n");

  174. 		printf("0: The droplet surface recedes and the droplet losses mass\n");

  175. 		printf("1: The droplet surface does not move and the droplet mass is constant\n");

  176. 		return(NULL);

  177. 	}

  178. 

  179. 	ier=parseNumber<double>(input, "dPdt", MAXBUFLEN, &data->dPdt);

  180. 

  181. 	ier=parseNumber<double>(input, "Rg", MAXBUFLEN, &data->Rg);

  182. 	if(data->Rg < 0.0){

  183. 		printf("Rg must be greater than 0");

  184. 		return(NULL);

  185. 	}

  186. 

  187. 	ier=parseNumber<int>   (input, "reflectProblem"	, MAXBUFLEN, &data->reflectProblem);

  188. 	if(data->reflectProblem!=0 && data->reflectProblem!=1){

  189. 		printf("Invalid entry for reflectProblem! Can be only 1 or 0.\n");

  190. 		return(NULL);

  191. 	}

  192. 	

  193. 	ier=parseNumber<double>(input, "mdot"		, MAXBUFLEN, &data->mdot);

  194. 	

  195. 	ier=parseNumber<double>(input, "initialTemperature", MAXBUFLEN,

  196. 			&data->initialTemperature);

  197. 	if(ier==-1 || data->initialTemperature<=0.0e0){

  198. 		printf("Enter positive initialTemperature in K!\n");

  199. 		return(NULL);

  200. 	}

  201. 

  202. 

  203. 	ier=parseNumber<double>(input, "initialPressure", MAXBUFLEN,

  204. 			&data->initialPressure);

  205. 	if(ier==-1 || data->initialTemperature<=0.0e0){

  206. 		printf("Enter positive initialPressure in atm!\n");

  207. 		return(NULL);

  208. 	}

  209. 

  210. 	ier=parseNumber<int>   (input, "metric"	, MAXBUFLEN, &data->metric);

  211. 	if(data->metric!=0 && data->metric!=1 && data->metric!=2){

  212. 		printf("Invalid entry for metric!\n");

  213. 		printf("0: Cartesian\n");

  214. 		printf("1: Cylindrical\n");

  215. 		printf("2: Spherical\n");

  216. 		return(NULL);

  217. 	}

  218. 

  219. 	ier=parseNumber<double>(input, "QDot", MAXBUFLEN, &data->maxQDot);

  220. 	if(ier==-1 && data->problemType==0){

  221. 		printf("Need to specify QDot for problemType 0!\n");

  222. 		return(NULL);

  223. 	}

  224. 

  225. 	ier=parseNumber<double>(input, "kernelSize", MAXBUFLEN, &data->kernelSize);

  226. 	if(ier==-1 && data->problemType==0){

  227. 		printf("Need to specify kernelSize for problemType 0!\n");

  228. 		return(NULL);

  229. 	}

  230. 

  231. 	ier=parseNumber<double>(input, "ignTime", MAXBUFLEN, &data->ignTime);

  232. 	if(ier==-1 && data->problemType==0){

  233. 		printf("Need to specify ignTime for problemType 0!\n");

  234. 		return(NULL);

  235. 	}

  236. 

  237. 	ier=parseNumber<double>(input, "mixingWidth", MAXBUFLEN,

  238. 			&data->mixingWidth);

  239. 	if(ier==-1){

  240. 		printf("Need to specify mixingWidth!\n");

  241. 		return(NULL);

  242. 	}

  243. 	

  244. 	ier=parseNumber<double>(input, "shift", MAXBUFLEN, &data->shift);

  245. 

  246. 	ier=parseNumber<double>(input, "firstRadius", MAXBUFLEN, &data->firstRadius);

  247. 	

  248. 	ier=parseNumber<double>(input, "wallTemperature", MAXBUFLEN, &data->wallTemperature);

  249. 

  250. 	ier=parseNumber<int>   (input, "dirichletInner"	, MAXBUFLEN,

  251. 			&data->dirichletInner);

  252. 	if(data->dirichletInner!=0 && data->dirichletInner!=1){

  253. 		printf("dirichletInner can either be 0 or 1!\n");

  254. 		return(NULL);

  255. 	}

  256. 

  257. 	ier=parseNumber<int>   (input, "dirichletOuter"	, MAXBUFLEN,

  258. 			&data->dirichletOuter);

  259. 	if(data->dirichletOuter!=0 && data->dirichletOuter!=1){

  260. 		printf("dirichletOuter can either be 0 or 1!\n");

  261. 		return(NULL);

  262. 	}

  263. 

  264. 	ier=parseNumber<int>   (input, "adaptiveGrid"	, MAXBUFLEN,

  265. 			&data->adaptiveGrid);

  266. 	if(ier==-1 || (data->adaptiveGrid!=0 && data->adaptiveGrid!=1)){

  267. 		printf("specify adaptiveGrid as 0 or 1!\n");

  268. 		return(NULL);

  269. 	}

  270. 

  271. 	ier=parseNumber<int>   (input, "moveGrid"	, MAXBUFLEN,

  272. 			&data->moveGrid);

  273. 	if(ier==-1 || (data->moveGrid!=0 && data->moveGrid!=1)){

  274. 		printf("specify moveGrid as 0 or 1!\n");

  275. 		return(NULL);

  276. 	}

  277. 

  278. 	ier=parseNumber<double>   (input, "isotherm"	, MAXBUFLEN,

  279. 			&data->isotherm);

  280. 	if(ier==-1){

  281. 		printf("specify temperature of isotherm!\n");

  282. 		return(NULL);

  283. 	}

  284. 

  285. 	ier=parseNumber<double>(input, "gridOffset", MAXBUFLEN, &data->gridOffset);

  286. 

  287. 	ier=parseNumber<int>   (input, "nSaves"	, MAXBUFLEN, &data->nSaves);

  288. 	if(data->nSaves<0 ){

  289. 		printf("nSaves must be greater than 0!\n");

  290. 		return(NULL);

  291. 	}

  292. 	

  293. 	ier=parseNumber<int>   (input, "writeRates"	, MAXBUFLEN,

  294. 			&data->writeRates);

  295. 	if(data->writeRates!=0 && data->writeRates!=1){

  296. 		printf("writeRates must either be 0 or 1!\n");

  297. 		return(NULL);

  298. 	}

  299. 

  300. 	ier=parseNumber<int>   (input, "writeEveryRegrid", MAXBUFLEN,

  301. 			&data->writeEveryRegrid);

  302. 	if(data->writeEveryRegrid!=0 && data->writeEveryRegrid!=1){

  303. 		printf("writeEveryRegrid must either be 0 or 1!\n");

  304. 		return(NULL);

  305. 	}

  306. 

  307. 	ier=parseNumber<int>   (input, "setConstraints"	, MAXBUFLEN,

  308. 			&data->setConstraints);

  309. 	if(data->setConstraints!=0 && data->setConstraints!=1){

  310. 		printf("setConstraints must either be 0 or 1!\n");

  311. 		return(NULL);

  312. 	}

  313. 

  314. 	ier=parseNumber<int>   (input, "suppressAlg"	, MAXBUFLEN,

  315. 			&data->suppressAlg);

  316. 	if(data->suppressAlg!=0 && data->suppressAlg!=1){

  317. 		printf("suppressAlg must either be 0 or 1!\n");

  318. 		return(NULL);

  319. 	}

  320. 

  321. 	ier=parseNumber<int>   (input, "dryRun"	, MAXBUFLEN,

  322. 			&data->dryRun);

  323. 	if(data->dryRun!=0 && data->dryRun!=1){

  324. 		printf("dryRun must either be 0 or 1!\n");

  325. 		return(NULL);

  326. 	}

  327. 

  328. 	ier=parseNumber<double>   (input, "finalTime"	, MAXBUFLEN,

  329. 			&data->finalTime);

  330. 

  331. 	ier=parseNumber<double>   (input, "relativeTolerance"	, MAXBUFLEN,

  332. 			&data->relativeTolerance);

  333. 	ier=parseNumber<double>   (input, "radiusTolerance"	, MAXBUFLEN,

  334. 			&data->radiusTolerance);

  335. 	ier=parseNumber<double>   (input, "temperatureTolerance", MAXBUFLEN,

  336. 			&data->temperatureTolerance);

  337. 	ier=parseNumber<double>   (input, "pressureTolerance", MAXBUFLEN,

  338. 			&data->pressureTolerance);

  339. 	ier=parseNumber<double>   (input, "massFractionTolerance", MAXBUFLEN,

  340. 			&data->massFractionTolerance);

  341. 	ier=parseNumber<double>   (input, "bathGasTolerance", MAXBUFLEN,

  342. 			&data->bathGasTolerance);

  343. 	ier=parseNumber<double>   (input, "MdotTolerance", MAXBUFLEN,

  344. 			&data->MdotTolerance);

  345. 

  346. 	char chem[MAXBUFLEN],mix[MAXBUFLEN],tran[MAXBUFLEN];

  347. 

  348. 	ier=parseNumber<char>(input, "chemistryFile"	, MAXBUFLEN, chem);

  349. 	if(ier==-1){

  350. 		printf("Enter chemistryFile!\n");

  351. 		return(NULL);

  352. 	}else{

  353. 		try{

  354. 			data->gas = new Cantera::IdealGasMix(chem);

  355. 			data->nsp=data->gas->nSpecies(); //assign no: of species

  356. 

  357.     		} catch (Cantera::CanteraError& err) {

  358. 		    printf("Error:\n");

  359. 		    printf("%s\n",err.what());

  360. 		    return(NULL);

  361.     		}

  362. 	}

  363. 

  364. //	/*No. of liquid phase species is determined by droplet type, either 1 or 2*/

  365. //	if(data->dropType==0){

  366. //		data->l_nsp=1;

  367. //	}else{

  368. //		data->l_nsp=2;

  369. //	}

  370. 

  371. 	ier=parseNumber<char>(input, "transportModel", MAXBUFLEN, tran);

  372. 	if(ier==-1){

  373. 		printf("Enter transportModel!\n");

  374. 		return(NULL);

  375. 	}else{

  376. 		try{

  377.   			data->trmix = Cantera::newTransportMgr(tran, data->gas);

  378. 		}catch (Cantera::CanteraError& err) {

  379. 			printf("Error:\n");

  380. 			printf("%s\n",err.what());

  381. 			return(NULL);

  382.     		}

  383. 	}

  384. 

  385. 	ier=parseNumber<char>(input, "mixtureComposition", MAXBUFLEN, mix);

  386. 	if(ier==-1){

  387. 		printf("Enter mixtureComposition!\n");

  388. 		return(NULL);

  389. 	}else{

  390. 		if(data->gas!=NULL){

  391. 			try{

  392. 				data->gas->setState_TPX(data->initialTemperature,

  393. 							data->initialPressure*Cantera::OneAtm,

  394. 							mix);

  395. 			}catch (Cantera::CanteraError& err) {

  396. 		    		printf("Error:\n");

  397. 		    		printf("%s\n",err.what());

  398. 		    		return(NULL);

  399.     			}

  400. 		}

  401. 	}

  402. 	

  403. 	//ier=parseNumber<char>(input, "dropletComposition", MAXBUFLEN, data->dropSpec);

  404. 	//if(ier==-1){

  405. 	//	printf("Enter composition of droplet!\n");

  406. 	//	return(NULL);

  407. 	//}

  408. 

  409.     //ier=parseDrop(input,"dropletComposition",data->dropSpec,data->dropMole,MAXBUFLEN);

  410.     //ier=parseDrop(input,"dropletDensity",data->dropSpec,data->dropDens,MAXBUFLEN);

  411. 

  412. 

  413. 	ier=parseNumber<int>   (input, "nThreads", MAXBUFLEN, &data->nThreads);

  414. 	if(data->nThreads<0 ){

  415. 		printf("nThreads must be greater than 0!\n");

  416. 		return(NULL);

  417. 	}

  418.     

  419. 

  420. 	ier=parseNumber<double>(input, "PCAD", MAXBUFLEN, &data->PCAD);

  421. 	// if(ier==-1){

  422. 	// 	printf("Enter PCAD value!\n");

  423. 	// 	return(NULL);

  424. 	// }

  425.     ier=parseNumber<double>(input,"RGTC", MAXBUFLEN, &data->RGTC);

  426. 	// if(ier==-1){

  427. 	// 	printf("Enter PGTC value!\n");

  428. 	// 	return(NULL);

  429. 	// }

  430.     ier=parseNumber<int>(input,"JJRG", MAXBUFLEN, &data->JJRG);

  431. 	// if(ier==-1){

  432. 	// 	printf("Enter JJRG value!\n");

  433. 	// 	return(NULL);

  434. 	// }

  435.     ier=parseNumber<double>(input,"deltaT", MAXBUFLEN, &data->deltaT);

  436. 	// if(ier==-1){

  437. 	// 	printf("Enter deltaT value!\n");

  438. 	// 	return(NULL);

  439. 	// }

  440. 	

  441.     data->nr=0;

  442. 	//data->np=data->nr+1;

  443. 	data->nt=data->nr+1;

  444. 	data->ny=data->nt+1;

  445. 	data->np=data->ny+data->nsp;

  446. 	data->nm=data->np+1;			

  447. 	data->nvar=data->nsp+4;		 //assign no: of variables (R,T,P,Mdot,nsp species)

  448. 

  449. 	// data->l_nr=0;

  450. 	// data->np=data->nr+1;

  451. 	// data->l_nt=data->l_nr+1;

  452. 	// data->l_ny=data->l_nt+1;

  453. 	// data->l_np=data->l_ny+data->l_nsp; //l_nsp is either 1 or 2

  454. 	// data->l_nm=data->l_np+1;			

  455. 	// data->l_nvar=data->l_nsp+4;

  456. 

  457. 	data->npts = data->l_npts+data->g_npts;

  458.     //std::vector<double> MW(data->dropMole.size());

  459.     for(size_t ii=0;ii< data->dropMole.size();ii++){

  460. //        printf("Print inside the loop.\n" ) ;

  461. //        printf("Print dropSpec component :%s! \n",data->dropSpec[0].c_str());

  462. //        size_t index = 0;

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

  464.             printf("spcies name :%s! \t",data->gas->speciesName(k-1).c_str());

  465.             if(data->gas->speciesName(k - 1)==data->dropSpec[ii] ) {

  466. //                index = k;

  467.                 double weight = data->gas->molecularWeight(k-1);

  468.                 printf("Molecular weight is: %.3f. \n",weight) ;

  469.                 data->MW.push_back(weight);

  470.             }

  471.         }

  472.     }

  473. //    printf("MW array size: %zu\n",data->MW.size());

  474.     for(auto & arg :data->MW){

  475.         printf("Molecular weight : %.3f. \n",arg) ;

  476.     }

  477.     // double massSum = 0.0; 

  478.     // for(int i=0;i<=1;i++){

  479.     //     massSum = massSum + data->dropMole[i] * MW[i];

  480.     // }

  481.     // data->dropRho = 0.0;

  482.     // for(int i=0;i<=1;i++){

  483.     //     data->dropMassFrac[i] = data->dropMole[i]*MW[i]/massSum;

  484.     //     data->dropRho = data->dropRho + data->dropMassFrac[i]*data->dropDens[i];

  485.     // }

  486. 	

  487.     // printf("Density of droplet is: %.3f [kg/m^3]\n",data->dropRho);

  488. 	// Mass of droplet

  489. 	//data->massDrop=1.0/3.0*pow(data->Rd,3)*997.0; //TODO: The density of the droplet should be a user input

  490.     //data->massDrop=1.0/3.0*pow(data->Rd,3)*684.00; //TODO: The density of the droplet(n-heptane) should be a user input

  491.     // data->massDrop=1.0/3.0*pow(data->Rd,3)*data->dropRho;

  492. 

  493. 	/* Update the gamma using UNIFAC method for propane and n-heptane */

  494. 	/*Note: gamma only need to be calculated when droplet type = 1, i.e. bi-component*/

  495. 	getGamma(data->dropMole,data->gamma); 

  496. //	printf("Gamma of Propane and n-heptane are %.3f and %.3f .\n",data->gamma[0],data->gamma[1]) ;

  497.     if(data->dropType == 0){

  498.         printf("Single component droplet problem, activity coeffcient is 1. \n");

  499.     }else{

  500.         for(auto & arg : data->gamma) {

  501.             printf("Initial activity coeffcient %.3f .\n",arg);

  502.         }

  503.     }

  504. 

  505. 	if(!data->adaptiveGrid){

  506. 		data->uniformGrid = new double [data->npts];

  507. 		// data->g_neq=data->g_nvar*data->g_npts; 

  508. 		// data->l_neq=data->l_nvar*data->l_npts;

  509. 		data->neq=data->npts*data->nvar;

  510. 	}

  511. 	else{

  512. 		data->grid=(UserGrid) malloc(sizeof *data->grid);

  513. 		ier=getGridSettings(input,data->grid);

  514. 		if(ier==-1)return(NULL);

  515. 

  516. 		ier=initializeGrid(data->grid);

  517. 		if(ier==-1)return(NULL);

  518. 

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

  520. 		if(ier==-1)return(NULL);

  521. 

  522. 		data->npts=data->grid->nPts;

  523. 		data->neq=data->nvar*data->npts;

  524. 	}

  525. 

  526. 	data->output=fopen("output.dat","w");

  527. 	data->globalOutput=fopen("globalOutput.dat","w");

  528. 	data->gridOutput=fopen("grid.dat","w");

  529. 	// data->timescaleOutput=fopen("timeScale.dat","w") ;

  530.     // data->rxnROPOutput=fopen("rxnROP.dat","w");

  531.     // data->spROPOutput=fopen("spROP.dat","w");

  532. 	//data->ratesOutput=fopen("rates.dat","w");

  533. 	

  534. 	data->innerMassFractions = new double [data->nsp];

  535. 	//data->wdot_mole = new double [data->nsp] ; 

  536. 	//data->wdot_mass = new double [data->nsp] ; 

  537. 	data->time_scale = new double [(data->npts) * (data->nsp)] ; 

  538. 	//data->MW = new double [data->nsp] ; 

  539. 

  540. 	return(data);

  541. }

  542. 

  543. void setSaneDefaults(UserData data){

  544. 	data->domainLength=1.0e-02;

  545. 	data->Rd=100.0e-06;

  546. 	data->constantPressure=1;

  547. 	data->problemType=1;

  548. 	data->quasiSteady=1;

  549. 	data->dPdt=0.0e0;

  550. 	data->reflectProblem=0;

  551. 	data->mdot=0.0;

  552. 	data->initialTemperature=300.0;

  553. 	data->initialPressure=1.0;

  554. 	data->metric=0;

  555. 	data->ignTime=1e-09;

  556. 	data->maxQDot=0.0e0;

  557. 	data->maxTemperature=300.0e0;

  558. 	data->mixingWidth=1e-04;

  559. 	data->shift=3.0e-03;

  560. 	data->firstRadius=1e-04;

  561. 	data->wallTemperature=298.0e0;

  562. 	data->dirichletInner=0;

  563. 	data->dirichletOuter=0;

  564. 	data->adaptiveGrid=0;

  565. 	data->moveGrid=0;

  566. 	data->gridOffset=0.0e0;

  567. 	data->Rg=1.0;

  568. 	data->isotherm=1000.0;

  569. 	data->nSaves=30;

  570. 	data->writeRates=0;

  571. 	data->writeEveryRegrid=0;

  572. 	data->relativeTolerance=1e-06;

  573. 	data->radiusTolerance=1e-08;

  574. 	data->temperatureTolerance=1e-06;

  575. 	data->pressureTolerance=1e-06;

  576. 	data->massFractionTolerance=1e-09;

  577. 	data->bathGasTolerance=1e-06;

  578. 	data->MdotTolerance=1e-09;

  579. 	data->finalTime=1e-02;

  580. 	data->tNow=0.0e0;

  581. 	data->setConstraints=0;

  582. 	data->suppressAlg=1;

  583. 	data->regrid=0;

  584. 	data->gridDirection=1;

  585. 	data->dryRun=0;

  586. 	data->nThreads=1;

  587. 

  588. 	data->flamePosition[0]=0.0e0;

  589. 	data->flamePosition[1]=0.0e0;

  590. 	data->flameTime[0]=0.0e0;

  591. 	data->flameTime[1]=0.0e0;

  592. 	data->nTimeSteps=0;

  593.     data->PCAD=0.75;

  594.     data->RGTC=1.0;

  595.     data->JJRG=0;

  596.     data->deltaT=200.0;

  597. }

  598. 

  599. 

  600. void getGamma(const double mole[],double gamma[]){

  601.     /* define relevant matrix */

  602.     Eigen::Matrix2d nu_ ;

  603.     nu_ << 2,1,2,5;

  604.     /* define relevant vectors */

  605.     Eigen::Vector2d R_(0.9011,0.6744),Q_(0.8480,0.5400);

  606.     Eigen::Vector2d r_,q_,x_,l_,ones_(1.0,1.0),v_,ksi_,gammaC_;

  607.     double sum_q, sum_r,sum_l;

  608. 

  609.     r_ = nu_ * R_ ;

  610.     q_ = nu_ * Q_ ;

  611.     x_ << mole[0], mole[1] ;

  612. 

  613.     sum_q = x_.dot(q_) ;

  614.     sum_r = x_.dot(r_) ;

  615. 

  616.     for(int i=0; i<v_.size() ; i++){

  617.         v_[i] = q_[i] * x_[i] / sum_q;

  618.         ksi_[i] = r_[i] * x_[i] / sum_r ;

  619.     }

  620. 

  621.     l_ = 5 * (r_ - q_) - (r_ - ones_) ;

  622.     sum_l = l_.dot(x_) ;

  623.     /* calculate the gamma_c  */

  624.     for(int i=0; i<v_.size() ; i++){

  625.         gammaC_[i] =std::exp( std::log(ksi_[i]/x_[i]) + 5*q_[i]*std::log(v_[i]/ksi_[i]) + l_[i] - ksi_[i]/x_[i]*sum_l );

  626.     }

  627.     /******* gamma_R for both propane and n-heptane are 0 **********/

  628.     /* return the gamma */

  629.     for(int i=0; i<v_.size(); i++){

  630.         gamma[i] = gammaC_[i] ;

  631.     }

  632. }

  633. 

  634. 

  635. void getGamma(const std::vector<double>& mole,std::vector<double>& gamma){

  636. 	size_t sz = mole.size();

  637. 	if (sz == 1){

  638. 		gamma.push_back(1); 

  639. 	}else{

  640.         gamma.reserve(mole.size());

  641. 		/* define relevant matrix */

  642. 		Eigen::Matrix2d nu_ ;

  643. 		nu_ << 2,1,2,5;

  644. 		/* define relevant vectors */

  645. 		Eigen::Vector2d R_(0.9011,0.6744),Q_(0.8480,0.5400);

  646. 		Eigen::Vector2d r_,q_,x_,l_,ones_(1.0,1.0),v_,ksi_,gammaC_;

  647. 		double sum_q, sum_r,sum_l;

  648. 

  649. 		r_ = nu_ * R_ ;

  650. 		q_ = nu_ * Q_ ;

  651. 		x_ << mole[0], mole[1] ;

  652. 

  653. 		sum_q = x_.dot(q_) ;

  654. 		sum_r = x_.dot(r_) ;

  655. 

  656. 		for(size_t i=0; i<v_.size() ; i++){

  657. 			v_[i] = q_[i] * x_[i] / sum_q;

  658. 			ksi_[i] = r_[i] * x_[i] / sum_r ;

  659. 		}

  660. 

  661. 		l_ = 5 * (r_ - q_) - (r_ - ones_) ;

  662. 		sum_l = l_.dot(x_) ;

  663. 		/* calculate the gamma_c  */

  664. 		for(size_t i=0; i<v_.size() ; i++){

  665. 			gammaC_[i] =std::exp( std::log(ksi_[i]/x_[i]) + 5*q_[i]*std::log(v_[i]/ksi_[i]) + l_[i] - ksi_[i]/x_[i]*sum_l );

  666. 		}

  667. 		/******* gamma_R for both propane and n-heptane are 0 **********/

  668. 		/* return the gamma */

  669. 		for(size_t i=0; i<v_.size(); i++){

  670. 			gamma.push_back(gammaC_[i]);

  671. 		}

  672. 	}

  673. 	

  674. }