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

4 anos atrás
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338
  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. }