Droplet Lagrangian Transient One-dimensional Reacting Code Implementation of both liquid and gas phase governing equations.
Nevar pievienot vairāk kā 25 tēmas Tēmai ir jāsākas ar burtu vai ciparu, tā var saturēt domu zīmes ('-') un var būt līdz 35 simboliem gara.

294 rindas
6.2KB

  1. #include "gridRoutines.h"
  2. #include <stdio.h>
  3. inline double l(const double* x,
  4. const double* a,
  5. const double* w,
  6. const double* fac,
  7. const int* refineLeft){
  8. if(*refineLeft==0){
  9. return(tanh(-*a*(*x+*w*100.0e0)));
  10. }else{
  11. double l ;
  12. l = tanh(-*a*(*x-*w*(*fac))) ;
  13. if(l>=0){
  14. return l*10.0;
  15. }else{
  16. //return(tanh(-*a*(*x-*w*(*fac))));
  17. return l;
  18. }
  19. }
  20. }
  21. //inline double l(const double* x,
  22. // const double* a,
  23. // const double* c,
  24. // const double* w,
  25. // const double* fac,
  26. // const int* refineLeft){
  27. // if(*refineLeft==0){
  28. // return(tanh(-*a*(*x+*w*100.0e0)));
  29. // }else{
  30. // double l ;
  31. // //l = tanh(-*a*(*x-*w*(*fac))) ;
  32. // l = tanh(-*a*(*x-*c));
  33. //
  34. // if(l>=0){
  35. // return l*10.0;
  36. // }else{
  37. // //return(tanh(-*a*(*x-*w*(*fac))));
  38. // return l;
  39. // }
  40. // }
  41. //}
  42. inline double r(const double* x,
  43. const double* a,
  44. const double* w,
  45. const double* fac,
  46. const int* refineRight){
  47. if(*refineRight==0){
  48. return(tanh(*a*(*x-(1.0e0+*w*100.0e0))));
  49. }else{
  50. return(tanh(*a*(*x-(1.0e0-*w*(*fac)))));
  51. }
  52. }
  53. inline double f(const double* x,
  54. const double* a,
  55. const double* c,
  56. const double* w){
  57. return(tanh(-*a*(*x-(*c+*w)))
  58. +tanh(-*a*((*x-1.0e0)-(*c+*w)))
  59. +tanh(-*a*((*x+1.0e0)-(*c+*w))));
  60. }
  61. inline double g(const double* x,
  62. const double* a,
  63. const double* c,
  64. const double* w){
  65. return(tanh(*a*(*x-(*c-*w)))
  66. +tanh(*a*((*x-1.0e0)-(*c-*w)))
  67. +tanh(*a*((*x+1.0e0)-(*c-*w))));
  68. }
  69. inline double rho(const double* x,
  70. const double* a,
  71. const double* c,
  72. const double* w,
  73. const double* mag,
  74. const double* leftFac,
  75. const double* rightFac,
  76. const int* refineLeft,
  77. const int* refineRight){
  78. return(((2.0e0+f(x,a,c,w)
  79. +g(x,a,c,w)
  80. +l(x,a,w,leftFac,refineLeft)
  81. // +l(x,a,c,w,leftFac,refineLeft)
  82. +r(x,a,w,rightFac,refineRight))*0.5e0)
  83. *(*mag-1.0e0)+1.0e0);
  84. }
  85. size_t maxPoints(const size_t basePts,
  86. const double* a,
  87. const double* w,
  88. const double* mag,
  89. const double* leftFac,
  90. const double* rightFac,
  91. const int* refineLeft,
  92. const int* refineRight){
  93. double dx=1.0e0/((double)(basePts)-1.0e0);
  94. double y=0.0e0;
  95. size_t i=0;
  96. double r=0.0e0;
  97. double t=0.5e0;
  98. while(y<=1.0e0){
  99. r=rho(&y,a,&t,w,mag,leftFac,rightFac,refineLeft,refineRight);
  100. y=y+(dx/r);
  101. i++;
  102. }
  103. return(i);
  104. }
  105. void fillGrid(const size_t* basePts,
  106. const size_t* nPts,
  107. const double* a,
  108. const double* c,
  109. const double* w,
  110. const double* mag,
  111. const double* leftFac,
  112. const double* rightFac,
  113. const int* refineLeft,
  114. const int* refineRight,
  115. double x[]){
  116. FILE* out;out=fopen("tmp.dat","w");
  117. double y=0.0e0;
  118. double r=0.0e0;
  119. double dx=1.0e0/((double)(*basePts)-1.0e0);
  120. for(size_t j=0;j<*nPts;j++){
  121. r=rho(&y,a,c,w,mag,leftFac,rightFac,refineLeft,refineRight); // Point density?
  122. fprintf(out, "%15.15e\n",dx/r); // writing number of points per section to tmp.dat
  123. y=y+(dx/r); // y is the total number of points?
  124. }
  125. fclose(out);
  126. double dxp[*nPts-1];
  127. for (size_t j = 0; j < *nPts; j++) {
  128. dxp[j]=0.0e0;
  129. }
  130. FILE* tmp;tmp=fopen("tmp.dat","r");
  131. char buf[MAXBUFLEN];
  132. size_t i=0;
  133. while (fgets(buf,MAXBUFLEN, tmp)!=NULL){
  134. sscanf(buf, "%lf", &y);
  135. dxp[i]=y;
  136. i++;
  137. }
  138. fclose(tmp);
  139. double sum=0.0e0;
  140. double err=0.0e0;
  141. double fix=0.0e0;
  142. double arr[*nPts-1] ;
  143. size_t halfboundII = 0;
  144. for(size_t j=0;j<*nPts-1;j++){
  145. sum+=dxp[j];
  146. arr[j]=sum;
  147. }
  148. for(size_t j=0;j<*nPts;j++){
  149. if(arr[j] > 0.5e0){
  150. halfboundII = j;
  151. break;
  152. }
  153. }
  154. err=1.0e0-sum;
  155. printf("sum before correction: %15.6e\n",sum);
  156. printf("err before correction: %15.6e\n",err);
  157. //fix=err/((double)(*nPts));
  158. fix = err/((double)(*nPts-1-(halfboundII+1)));
  159. sum=0.0e0;
  160. for(size_t j=0;j<*nPts-1;j++){
  161. // dxp[j]+=fix;
  162. if(j>halfboundII){
  163. dxp[j]=dxp[j]+fix ;
  164. }
  165. sum+=dxp[j];
  166. }
  167. err=1.0e0-sum;
  168. printf("sum after correction:%15.6e\n",sum);
  169. printf("err after correction: %15.6e\n",err);
  170. x[0]=0.0e0;
  171. for(size_t j=0;j<*nPts-1;j++){
  172. x[j+1]=x[j]+dxp[j];
  173. }
  174. x[*nPts-1]=1.0e0;
  175. }
  176. double safePosition(double c, double w){
  177. if(c<w){
  178. return(w);
  179. }
  180. else if(c>1.0e0-w){
  181. return(1.0e0-w);
  182. }
  183. else{
  184. return(c);
  185. }
  186. }
  187. int reGrid(UserGrid grid, double position){
  188. printf("before regrid: %ld\n", grid->nPts);
  189. double xx[grid->nPts];
  190. fillGrid(&grid->basePts,
  191. &grid->nPts,
  192. &grid->a,
  193. &position,
  194. &grid->w,
  195. &grid->mag,
  196. &grid->leftFac,
  197. &grid->rightFac,
  198. &grid->refineLeft,
  199. &grid->refineRight,
  200. xx);
  201. for (size_t i = 0; i < grid->nPts; i++) {
  202. grid->x[i]=xx[i];
  203. }
  204. return(0);
  205. }
  206. void storeGrid(const double* x, double *y, const size_t nPts){
  207. for(size_t i=0;i<nPts;i++){
  208. y[i]=x[i];
  209. }
  210. }
  211. int initializeGrid(UserGrid grid){
  212. grid->nPts=maxPoints(grid->basePts,
  213. &grid->a,
  214. &grid->w,
  215. &grid->mag,
  216. &grid->leftFac,
  217. &grid->rightFac,
  218. &grid->refineLeft,
  219. &grid->refineRight);
  220. printf("nPts: %ld\n",grid->nPts);
  221. grid->leastMove=grid->w;
  222. grid->x = new double [grid->nPts];
  223. grid->xOld = new double [grid->nPts];
  224. for (size_t i = 0; i < grid->nPts; i++) {
  225. grid->x[i]=0.0e0;
  226. grid->xOld[i]=0.0e0;
  227. }
  228. return(0);
  229. }
  230. int getGridSettings(FILE *input, UserGrid grid){
  231. int ier=0;
  232. ier=parseNumber<size_t>(input, "basePts" , MAXBUFLEN, &grid->basePts);
  233. if(ier==-1)return(-1);
  234. ier=parseNumber<double>(input, "gridDensitySlope", MAXBUFLEN, &grid->a);
  235. if(ier==-1)return(-1);
  236. ier=parseNumber<double>(input, "fineGridHalfWidth", MAXBUFLEN, &grid->w);
  237. if(ier==-1)return(-1);
  238. ier=parseNumber<double>(input, "gridRefinement", MAXBUFLEN, &grid->mag);
  239. if(ier==-1)return(-1);
  240. ier=parseNumber<double>(input, "leftRefineFactor", MAXBUFLEN, &grid->leftFac);
  241. if(ier==-1)return(-1);
  242. ier=parseNumber<double>(input, "rightRefineFactor", MAXBUFLEN, &grid->rightFac);
  243. if(ier==-1)return(-1);
  244. ier=parseNumber<int>(input, "refineLeft" , MAXBUFLEN, &grid->refineLeft);
  245. if(ier==-1)return(-1);
  246. ier=parseNumber<int>(input, "refineRight" , MAXBUFLEN, &grid->refineRight);
  247. if(ier==-1)return(-1);
  248. ier=parseNumber<double>(input, "position" , MAXBUFLEN, &grid->position);
  249. if(ier==-1)return(-1);
  250. return(0);
  251. }