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