A low Mach, 1D, reacting flow code.
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

245 lignes
5.1KB

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