1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscsys.h> 4 5 extern int wrap_revolve(int* check,int* capo,int* fine,int *snaps_in,int* info); 6 7 typedef struct _StackElement { 8 PetscInt stepnum; 9 Vec X; 10 Vec *Y; 11 PetscReal time; 12 PetscReal timeprev; 13 } *StackElement; 14 15 typedef struct _Stack { 16 PetscBool usecontroller; 17 PetscBool reverseonestep; 18 PetscInt stepsleft; 19 PetscInt check; 20 PetscInt capo; 21 PetscInt fine; 22 PetscInt info; 23 PetscInt top; /* The top of the stack */ 24 PetscInt maxelements; /* The maximum stack size */ 25 PetscInt numY; 26 MPI_Comm comm; 27 StackElement *stack; /* The storage */ 28 } Stack; 29 30 static PetscErrorCode StackCreate(MPI_Comm,Stack *,PetscInt,PetscInt); 31 static PetscErrorCode StackDestroy(Stack*); 32 static PetscErrorCode StackPush(Stack*,StackElement); 33 static PetscErrorCode StackPop(Stack*,StackElement*); 34 static PetscErrorCode StackTop(Stack*,StackElement*); 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "StackCreate" 38 static PetscErrorCode StackCreate(MPI_Comm comm,Stack *s,PetscInt size,PetscInt ny) 39 { 40 PetscErrorCode ierr; 41 42 PetscFunctionBegin; 43 s->top = -1; 44 s->maxelements = size; 45 s->comm = comm; 46 s->numY = ny; 47 48 ierr = PetscMalloc1(s->maxelements*sizeof(StackElement),&s->stack);CHKERRQ(ierr); 49 ierr = PetscMemzero(s->stack,s->maxelements*sizeof(StackElement));CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "StackDestroy" 55 static PetscErrorCode StackDestroy(Stack *s) 56 { 57 PetscInt i; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 if (s->top>-1) { 62 for (i=0;i<=s->top;i++) { 63 ierr = VecDestroy(&s->stack[i]->X);CHKERRQ(ierr); 64 ierr = VecDestroyVecs(s->numY,&s->stack[i]->Y);CHKERRQ(ierr); 65 ierr = PetscFree(s->stack[i]);CHKERRQ(ierr); 66 } 67 } 68 ierr = PetscFree(s->stack);CHKERRQ(ierr); 69 ierr = PetscFree(s);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "StackPush" 75 static PetscErrorCode StackPush(Stack *s,StackElement e) 76 { 77 PetscFunctionBegin; 78 if (s->top+1 >= s->maxelements) SETERRQ1(s->comm,PETSC_ERR_MEMC,"Maximum stack size (%D) exceeded",s->maxelements); 79 s->stack[++s->top] = e; 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "StackPop" 85 static PetscErrorCode StackPop(Stack *s,StackElement *e) 86 { 87 PetscFunctionBegin; 88 if (s->top == -1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Emptry stack"); 89 *e = s->stack[s->top--]; 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "StackTop" 95 static PetscErrorCode StackTop(Stack *s,StackElement *e) 96 { 97 PetscFunctionBegin; 98 *e = s->stack[s->top]; 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "TSTrajectorySet_Memory" 104 PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 105 { 106 PetscInt ns,i; 107 Vec *Y; 108 PetscReal timeprev; 109 StackElement e; 110 Stack *s = (Stack*)tj->data; 111 PetscInt whattodo,oldcapo; 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 if (stepnum<s->top) SETERRQ(s->comm,PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 116 117 if (s->usecontroller) { 118 printf("left %d\n",s->stepsleft); 119 if (s->reverseonestep) PetscFunctionReturn(0); 120 if (s->stepsleft==0) { /* let the controller determine what to do next */ 121 s->capo = stepnum; 122 oldcapo = s->capo; 123 printf("check=%d,capo=%d,fine=%d,snaps_in=%d\t",s->check,s->capo,s->fine,s->maxelements); 124 whattodo = wrap_revolve(&s->check,&s->capo,&s->fine,&s->maxelements,&s->info); 125 switch(whattodo) { 126 case 1: 127 printf("Advance from %d to %d.\n", oldcapo, s->capo); 128 break; 129 case 2: 130 printf("Store in checkpoint number %d\n",s->check); 131 break; 132 case 3: 133 printf("First turn: Initialize adjoints and reverse first step.\n"); 134 break; 135 case 4: 136 printf("Forward and reverse one step.\n"); 137 break; 138 case 5: 139 printf("Restore in checkpoint number %d\n",s->check); 140 break; 141 case -1: 142 printf("Error!"); 143 break; 144 } 145 if (whattodo==-1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Error in the controller"); 146 if (whattodo==1) { 147 s->stepsleft = s->capo-oldcapo-1; 148 PetscFunctionReturn(0); /* do not need to checkpoint */ 149 } 150 if (whattodo==3 || whattodo==4) { 151 s->reverseonestep = PETSC_TRUE; 152 PetscFunctionReturn(0); 153 } 154 if (whattodo==5) { 155 printf("check=%d,capo=%d,fine=%d,snaps_in=%d\t",s->check,s->capo,s->fine,s->maxelements); 156 oldcapo = s->capo; 157 whattodo = wrap_revolve(&s->check,&s->capo,&s->fine,&s->maxelements,&s->info); 158 if (whattodo!=1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Something weird happened"); 159 s->stepsleft = s->capo-oldcapo; 160 switch(whattodo) { 161 case 1: 162 printf("5Advance from %d to %d.\n", oldcapo, s->capo); 163 break; 164 case 2: 165 printf("5Store in checkpoint number %d\n",s->check); 166 break; 167 case 3: 168 printf("5First turn: Initialize adjoints and reverse first step.\n"); 169 break; 170 case 4: 171 printf("5Forward and reverse one step.\n"); 172 break; 173 case 5: 174 printf("5Restore in checkpoint number %d\n",s->check); 175 break; 176 case -1: 177 printf("Error!"); 178 break; 179 } 180 PetscFunctionReturn(0); 181 } 182 if (whattodo==2) { 183 printf("check=%d,capo=%d,fine=%d,snaps_in=%d\t",s->check,s->capo,s->fine,s->maxelements); 184 oldcapo = s->capo; 185 whattodo = wrap_revolve(&s->check,&s->capo,&s->fine,&s->maxelements,&s->info); 186 if (whattodo!=1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Something weird happened"); 187 if (stepnum==0) s->stepsleft = s->capo-oldcapo-1; 188 else s->stepsleft = s->capo-oldcapo-1; 189 switch(whattodo) { 190 case 1: 191 printf("2Advance from %d to %d.\n", oldcapo, s->capo); 192 break; 193 case 2: 194 printf("2Store in checkpoint number %d\n",s->check); 195 break; 196 case 3: 197 printf("2First turn: Initialize adjoints and reverse first step.\n"); 198 break; 199 case 4: 200 printf("2Forward and reverse one step.\n"); 201 break; 202 case 5: 203 printf("2Restore in checkpoint number %d\n",s->check); 204 break; 205 case -1: 206 printf("Error!"); 207 break; 208 } 209 } 210 } else { /* advance s->stepsleft time steps without checkpointing */ 211 s->stepsleft--; 212 PetscFunctionReturn(0); 213 } 214 } 215 216 /* checkpoint to memmory */ 217 if (stepnum==s->top) { /* overwrite the top checkpoint */ 218 ierr = StackTop(s,&e); 219 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 220 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 221 for (i=0;i<ns;i++) { 222 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 223 } 224 e->stepnum = stepnum; 225 e->time = time; 226 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 227 e->timeprev = timeprev; 228 } else { 229 ierr = PetscCalloc1(1,&e); 230 ierr = VecDuplicate(X,&e->X);CHKERRQ(ierr); 231 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 232 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 233 ierr = VecDuplicateVecs(Y[0],ns,&e->Y);CHKERRQ(ierr); 234 for (i=0;i<ns;i++) { 235 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 236 } 237 e->stepnum = stepnum; 238 e->time = time; 239 if (stepnum == 0) { 240 e->timeprev = e->time - ts->time_step; /* for consistency */ 241 } else { 242 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 243 e->timeprev = timeprev; 244 } 245 ierr = StackPush(s,e);CHKERRQ(ierr); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "TSTrajectoryGet_Memory" 252 PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 253 { 254 Vec *Y; 255 PetscInt nr,i; 256 StackElement e; 257 Stack *s = (Stack*)tj->data; 258 PetscReal stepsize; 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 if (s->usecontroller && s->reverseonestep) { 263 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 264 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); /* go backward */ 265 s->reverseonestep = PETSC_FALSE; 266 printf("reverse step %d\n",stepnum); 267 PetscFunctionReturn(0); 268 } 269 270 /* restore a checkpoint */ 271 ierr = StackTop(s,&e);CHKERRQ(ierr); 272 ierr = VecCopy(e->X,ts->vec_sol);CHKERRQ(ierr); 273 ierr = TSGetStages(ts,&nr,&Y);CHKERRQ(ierr); 274 for (i=0;i<nr ;i++) { 275 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 276 } 277 *t = e->time; 278 279 printf("read stepnum=%d to get %d\n",e->stepnum,stepnum); 280 if (e->stepnum < stepnum) { /* need recomputation */ 281 PetscInt whattodo; 282 s->capo = stepnum; 283 printf("check=%d,capo=%d,fine=%d,snaps_in=%d\t",s->check,s->capo,s->fine,s->maxelements); 284 whattodo = wrap_revolve(&s->check,&s->capo,&s->fine,&s->maxelements,&s->info); 285 switch(whattodo) { 286 case 5: 287 printf("R:Restore in checkpoint number %d\n",s->check); 288 break; 289 case -1: 290 printf("Error!"); 291 break; 292 } 293 ierr = TSSetTimeStep(ts,(*t)-e->timeprev);CHKERRQ(ierr); 294 /* reset ts context */ 295 PetscInt steps = ts->steps; 296 ts->steps = e->stepnum; 297 ts->ptime = e->time; 298 ts->ptime_prev = e->timeprev; 299 printf("need to recompute %d steps\n",stepnum-e->stepnum); 300 for (i=e->stepnum;i<stepnum;i++) { /* assume fixed step size */ 301 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 302 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 303 ierr = TSStep(ts);CHKERRQ(ierr); 304 if (ts->event) { 305 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 306 } 307 if (!ts->steprollback) { 308 ierr = TSPostStep(ts);CHKERRQ(ierr); 309 } 310 } 311 if (!s->reverseonestep) printf("error\n"); 312 ts->steps = steps; 313 ts->total_steps = stepnum; 314 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 315 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); /* go backward */ 316 if (stepnum-e->stepnum==1) { 317 ierr = StackPop(s,&e);CHKERRQ(ierr); 318 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 319 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 320 ierr = PetscFree(e);CHKERRQ(ierr); 321 } 322 s->reverseonestep = PETSC_FALSE; 323 printf("reverse step %d\n",stepnum); 324 } else if (e->stepnum == stepnum) { 325 ierr = TSSetTimeStep(ts,-(*t)+e->timeprev);CHKERRQ(ierr); /* go backward */ 326 ierr = StackPop(s,&e);CHKERRQ(ierr); 327 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 328 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 329 ierr = PetscFree(e);CHKERRQ(ierr); 330 } else { 331 SETERRQ2(s->comm,PETSC_ERR_ARG_OUTOFRANGE,"The current step no. is %D, but the step number at top of the stack is %D",stepnum,e->stepnum); 332 } 333 334 PetscFunctionReturn(0); 335 } 336 337 #undef __FUNCT__ 338 #define __FUNCT__ "TSTrajectoryDestroy_Memory" 339 PETSC_EXTERN PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 340 { 341 Stack *s = (Stack*)tj->data; 342 PetscErrorCode ierr; 343 344 PetscFunctionBegin; 345 ierr = StackDestroy(s);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 /*MC 350 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 351 352 Level: intermediate 353 354 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 355 356 M*/ 357 #undef __FUNCT__ 358 #define __FUNCT__ "TSTrajectoryCreate_Memory" 359 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 360 { 361 PetscInt nr,maxsteps; 362 Stack *s; 363 PetscErrorCode ierr; 364 365 PetscFunctionBegin; 366 tj->ops->set = TSTrajectorySet_Memory; 367 tj->ops->get = TSTrajectoryGet_Memory; 368 tj->ops->destroy = TSTrajectoryDestroy_Memory; 369 370 ierr = PetscCalloc1(1,&s); 371 s->maxelements = 3; /* will be provided by users */ 372 ierr = TSGetStages(ts,&nr,PETSC_IGNORE);CHKERRQ(ierr); 373 374 maxsteps = PetscMin(ts->max_steps,(PetscInt)(ceil(ts->max_time/ts->time_step))); 375 if (s->maxelements-1<maxsteps) { /* Need to use a controller */ 376 s->usecontroller = PETSC_TRUE; 377 s->reverseonestep = PETSC_FALSE; 378 s->check = -1; 379 s->capo = 0; 380 s->fine = maxsteps; 381 s->info = 2; 382 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,s->maxelements,nr);CHKERRQ(ierr); 383 } else { /* Enough space for checkpointing all time steps */ 384 s->usecontroller = PETSC_FALSE; 385 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,ts->max_steps+1,nr);CHKERRQ(ierr); 386 } 387 tj->data = s; 388 PetscFunctionReturn(0); 389 } 390