Add the HRR related function and parameters into LTORC
Du kannst nicht mehr als 25 Themen auswählen Themen müssen entweder mit einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

vor 4 Jahren
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. }