Add the HRR related function and parameters into LTORC
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.

339 lignes
9.2KB

  1. /*
  2. _____ ___ ____ ____
  3. |_ _/ _ \| _ \ / ___|
  4. | || | | | |_) | |
  5. | || |_| | _ <| |___
  6. |_| \___/|_| \_\\____|
  7. */
  8. #include "UserData.h"
  9. #include "solution.h"
  10. #include "residue.h"
  11. #include "macros.h"
  12. #include "timing.h"
  13. #include <ida/ida.h>
  14. #include <ida/ida_direct.h>
  15. #include <sunmatrix/sunmatrix_band.h>
  16. #include <sunlinsol/sunlinsol_lapackband.h>
  17. //#include <ida/ida_band.h>
  18. static int check_flag(void *flagvalue,
  19. const char *funcname,
  20. int opt);
  21. void freeAtLast(void* mem, N_Vector *y,
  22. N_Vector *ydot,
  23. N_Vector *res,
  24. N_Vector *id,
  25. N_Vector *atolv,
  26. N_Vector *constraints,UserData data);
  27. int main(){
  28. FILE *input;input=fopen("input.dat","r");
  29. UserData data;data=NULL;data=allocateUserData(input);
  30. fclose(input);
  31. data->clockStart=get_wall_time();
  32. if(data==NULL){
  33. printf("check input file!\n");
  34. freeUserData(data);
  35. return(-1);
  36. }
  37. long int ier,mu,ml,count,netf,ncfn,njevals,nrevals;
  38. realtype tNow,*atolvdata,*constraintsdata,finalTime,tolsfac;
  39. N_Vector y,ydot,id,res,atolv,constraints;
  40. y=ydot=id=res=atolv=constraints=NULL;
  41. ier=allocateSolution(data->neq,data->nThreads,&y,&ydot,&id,&res,&atolv,&constraints);
  42. ier=setAlgebraicVariables(&id,data);
  43. ier=setInitialCondition(&y,&ydot,data);
  44. if(ier==-1)return(-1);
  45. tNow=data->tNow;
  46. finalTime=data->finalTime;
  47. double* ydata;
  48. double* ydotdata;
  49. ydata = N_VGetArrayPointer_OpenMP(y);
  50. ydotdata = N_VGetArrayPointer_OpenMP(ydot);
  51. void *mem;mem=NULL;mem = IDACreate();
  52. ier = IDASetUserData(mem, data);
  53. ier = IDASetId(mem, id);
  54. ier = IDAInit(mem, residue, tNow, y, ydot);
  55. // Atol array
  56. atolvdata = N_VGetArrayPointer_OpenMP(atolv);
  57. for (size_t i = 1; i <=data->npts; i++) {
  58. atolT(i) = data->temperatureTolerance;
  59. for (size_t k = 1; k <=data->nsp; k++) {
  60. if(k!=data->k_bath){
  61. atolY(i,k) = data->massFractionTolerance;
  62. }
  63. else{
  64. atolY(i,k) = data->bathGasTolerance;
  65. }
  66. }
  67. atolR(i) = data->radiusTolerance;
  68. atolP(i) = data->pressureTolerance;
  69. }
  70. ier = IDASVtolerances(mem, data->relativeTolerance, atolv);
  71. mu = 2*data->nvar; ml = mu;
  72. SUNMatrix A; A=NULL;
  73. A=SUNBandMatrix(data->neq,mu,ml,mu+ml);
  74. SUNLinearSolver LS; LS=NULL;
  75. //LS=SUNBandLinearSolver(y,A);
  76. LS=SUNLapackBand(y,A);
  77. ier=IDADlsSetLinearSolver(mem,LS,A);
  78. //ier = IDABand(mem, data->neq, mu, ml);
  79. constraintsdata = N_VGetArrayPointer_OpenMP(constraints);
  80. if(data->setConstraints){
  81. for (size_t i = 1; i <=data->npts; i++) {
  82. for (size_t k = 1; k <=data->nsp; k++) {
  83. constraintsY(i,k) = ONE;
  84. }
  85. }
  86. ier=IDASetConstraints(mem, constraints);
  87. }
  88. ier = IDASetSuppressAlg(mem, data->suppressAlg);
  89. if(check_flag(&ier, "IDASetSuppressAlg", 1)) return(1);
  90. //ier= IDASetMaxNumStepsIC(mem, 1);
  91. //ier= IDASetMaxNumJacsIC(mem,8);
  92. //ier= IDASetMaxNumItersIC(mem,100);
  93. //ier= IDASetMaxBacksIC(mem,2000);
  94. //ier = IDASetLineSearchOffIC(mem,SUNTRUE);
  95. if(!data->dryRun){
  96. printf("Calculating Initial Conditions:\n");
  97. printf("Cross your fingers...\n");
  98. ier = IDACalcIC(mem, IDA_YA_YDP_INIT, 1e-05*finalTime);
  99. //If at first IDACalcIC doesn't succeed, try, try, try again:
  100. for (int i = 0; i < 10; i++) {
  101. ier = IDACalcIC(mem, IDA_YA_YDP_INIT, (1e-01+pow(10,i)*finalTime));
  102. if(ier==0){
  103. break;
  104. }
  105. }
  106. //...or else cry again :(
  107. if(check_flag(&ier, "IDACalcIC", 1)){
  108. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  109. return(-1);
  110. }else{
  111. printf("Initial (Consistent) Conditions Calculated!\n");
  112. }
  113. ier = IDASetInitStep(mem,1e-12);
  114. }
  115. printSpaceTimeHeader(data);
  116. printGlobalHeader(data);
  117. printSpaceTimeOutput(tNow, &y, data->output, data);
  118. printSpaceTimeOutput(tNow, &y, data->gridOutput, data);
  119. if(!data->dryRun){
  120. count=0;
  121. double dt=1e-08;
  122. double t1=0.0e0;
  123. double xOld=0.0e0;
  124. double x=0.0e0;
  125. double dx=0.0e0;
  126. double dxMin=1.0e0;
  127. double dxRatio=dx/dxMin;
  128. double writeTime=tNow;
  129. int move=0;
  130. int kcur=0;
  131. if(data->adaptiveGrid){
  132. dxMin=data->grid->leastMove;
  133. //xOld=maxCurvPosition(ydata, data->nt, data->nvar,
  134. // data->grid->x, data->npts);
  135. xOld=isothermPosition(ydata, data->isotherm, data->nt,
  136. data->nvar, data->grid->x, data->npts);
  137. }
  138. while (tNow<=finalTime) {
  139. t1=tNow;
  140. ier = IDASolve(mem, finalTime, &tNow, y, ydot, IDA_ONE_STEP);
  141. if(check_flag(&ier, "IDASolve", 1)){
  142. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  143. return(-1);
  144. }
  145. dt=tNow-t1;
  146. ier = IDAGetCurrentOrder(mem, &kcur);
  147. if(data->adaptiveGrid==1 && data->moveGrid==1){
  148. //x=maxCurvPosition(ydata, data->nt, data->nvar,
  149. // data->grid->x, data->npts);
  150. x=isothermPosition(ydata, data->isotherm, data->nt,
  151. data->nvar, data->grid->x, data->npts);
  152. //x=maxGradPosition(ydata, data->nt, data->nvar,
  153. // data->grid->x, data->npts);
  154. dx=x-xOld;
  155. if(dx*dxMin>0.0e0){
  156. move=1;
  157. }else{
  158. move=-1;
  159. }
  160. //if(fabs(dx)>=dxMin && x+(double)(move)*0.5e0*dxMin<=1.0e0){
  161. dxRatio=fabs(dx/dxMin);
  162. if(dxRatio>=1.0e0 && dxRatio<=2.0e0){
  163. printf("Regridding!\n");
  164. data->regrid=1;
  165. //printSpaceTimeOutput(tNow, &y, data->gridOutput, data);
  166. ier=reGrid(data->grid, x+(double)(move)*0.5e0*dxMin);
  167. if(ier==-1){
  168. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  169. return(-1);
  170. }
  171. updateSolution(ydata, ydotdata, data->nvar,
  172. data->grid->xOld,data->grid->x,data->npts);
  173. storeGrid(data->grid->x,data->grid->xOld,data->npts);
  174. xOld=x;
  175. printf("Regrid Complete! Restarting Problem at %15.6e s\n",tNow);
  176. ier = IDAReInit(mem, tNow, y, ydot);
  177. if(check_flag(&ier, "IDAReInit", 1)){
  178. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  179. return(-1);
  180. }
  181. ier = IDASetInitStep(mem,1e-01*dt);
  182. printf("Reinitialized! Calculating Initial Conditions:\n");
  183. printf("Cross your fingers...\n");
  184. ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e-01*dt);
  185. if(check_flag(&ier, "IDACalcIC", 1)){
  186. ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e+01*dt);
  187. }
  188. //Every once in a while, for reasons
  189. //that befuddle this mathematically
  190. //lesser author, IDACalcIC fails. Here,
  191. //I desperately try to make it converge
  192. //again by sampling increasingly larger
  193. //time-steps:
  194. for (int i = 0; i < 10; i++) {
  195. ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+(1e-01+pow(10,i)*dt));
  196. if(ier==0){
  197. break;
  198. }
  199. }
  200. //Failure :( Back to the drawing board:
  201. if(check_flag(&ier, "IDACalcIC", 1)){
  202. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  203. return(-1);
  204. }
  205. printf("Initial (Consistent) Conditions Calculated!\n");
  206. printf("------------------------------------------\n\n");
  207. //if(data->writeEveryRegrid){
  208. }
  209. }
  210. if(count%data->nSaves==0 && !data->writeEveryRegrid){
  211. printSpaceTimeOutput(tNow, &y, data->output, data);
  212. //if(data->writeRates){
  213. // printSpaceTimeRates(tNow, ydot, data);
  214. //}
  215. }
  216. /*Print global variables only if time-step is of high order!*/
  217. if(data->nTimeSteps==0){
  218. data->flamePosition[0]=0.0e0;
  219. data->flamePosition[1]=0.0e0;
  220. data->flameTime[0]=tNow;
  221. data->flameTime[1]=tNow;
  222. }
  223. if(kcur>=2 ){
  224. printGlobalVariables(tNow,&y,&ydot,data);
  225. }
  226. if(tNow>writeTime){
  227. writeTime+=data->writeDeltaT;
  228. printSpaceTimeOutput(tNow, &y, data->output, data);
  229. FILE* fp;
  230. fp=fopen("restart.bin","w");
  231. writeRestart(tNow, &y, &ydot, fp, data);
  232. fclose(fp);
  233. }
  234. ier = IDAGetNumErrTestFails(mem, &netf);
  235. ier = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
  236. ier = IDADlsGetNumJacEvals(mem, &njevals);
  237. ier = IDADlsGetNumResEvals(mem, &nrevals);
  238. printf("etf = %ld ,"
  239. "nlsf= %ld ,"
  240. "J=%ld ,"
  241. "R=%ld ,"
  242. "o=%d ,",netf, ncfn, njevals, nrevals, kcur);
  243. printf("Time=%15.6e s,",tNow);
  244. printf("frac: %15.6e\n",dxRatio);
  245. count++;
  246. data->nTimeSteps=count;
  247. }
  248. }
  249. SUNLinSolFree(LS);
  250. SUNMatDestroy(A);
  251. freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
  252. return(0);
  253. }
  254. void freeAtLast(void* mem,
  255. N_Vector *y,
  256. N_Vector *ydot,
  257. N_Vector *res,
  258. N_Vector *id,
  259. N_Vector *atolv,
  260. N_Vector *constraints,UserData data){
  261. IDAFree(&mem);
  262. freeSolution(y,ydot,res,id,atolv,constraints);
  263. freeUserData(data);
  264. }
  265. static int check_flag(void *flagvalue, const char *funcname, int opt)
  266. {
  267. int *errflag;
  268. /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
  269. if (opt == 0 && flagvalue == NULL) {
  270. fprintf(stderr,
  271. "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
  272. funcname);
  273. return(1);
  274. }
  275. /* Check if flag < 0 */
  276. else if (opt == 1) {
  277. errflag = (int *) flagvalue;
  278. if (*errflag < 0) {
  279. fprintf(stderr,
  280. "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
  281. funcname, *errflag);
  282. return(1);
  283. }
  284. }
  285. /* Check if function returned NULL pointer - no memory allocated */
  286. else if (opt == 2 && flagvalue == NULL) {
  287. fprintf(stderr,
  288. "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
  289. funcname);
  290. return(1);
  291. }
  292. return(0);
  293. }