Actual source code: traj.c

petsc-master 2020-01-28
Report Typos and Errors
  1:  #include <petsc/private/tsimpl.h>
  2:  #include <petsc/private/tshistoryimpl.h>
  3:  #include <petscdm.h>

  5: PetscFunctionList TSTrajectoryList              = NULL;
  6: PetscBool         TSTrajectoryRegisterAllCalled = PETSC_FALSE;
  7: PetscClassId      TSTRAJECTORY_CLASSID;
  8: PetscLogEvent     TSTrajectory_Set, TSTrajectory_Get, TSTrajectory_GetVecs;

 10: /*@C
 11:   TSTrajectoryRegister - Adds a way of storing trajectories to the TS package

 13:   Not Collective

 15:   Input Parameters:
 16: + name        - the name of a new user-defined creation routine
 17: - create_func - the creation routine itself

 19:   Notes:
 20:   TSTrajectoryRegister() may be called multiple times to add several user-defined tses.

 22:   Level: developer

 24: .seealso: TSTrajectoryRegisterAll()
 25: @*/
 26: PetscErrorCode TSTrajectoryRegister(const char sname[],PetscErrorCode (*function)(TSTrajectory,TS))
 27: {

 31:   PetscFunctionListAdd(&TSTrajectoryList,sname,function);
 32:   return(0);
 33: }

 35: /*@
 36:   TSTrajectorySet - Sets a vector of state in the trajectory object

 38:   Collective on TSTrajectory

 40:   Input Parameters:
 41: + tj      - the trajectory object
 42: . ts      - the time stepper object (optional)
 43: . stepnum - the step number
 44: . time    - the current time
 45: - X       - the current solution

 47:   Level: developer

 49:   Notes: Usually one does not call this routine, it is called automatically during TSSolve()

 51: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectoryGet(), TSTrajectoryGetVecs()
 52: @*/
 53: PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
 54: {

 58:   if (!tj) return(0);
 64:   if (!tj->ops->set) SETERRQ1(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"TSTrajectory type %s",((PetscObject)tj)->type_name);
 65:   if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
 66:   if (tj->monitor) {
 67:     PetscViewerASCIIPrintf(tj->monitor,"TSTrajectorySet: stepnum %D, time %g (stages %D)\n",stepnum,(double)time,(PetscInt)!tj->solution_only);
 68:   }
 69:   PetscLogEventBegin(TSTrajectory_Set,tj,ts,0,0);
 70:   (*tj->ops->set)(tj,ts,stepnum,time,X);
 71:   PetscLogEventEnd(TSTrajectory_Set,tj,ts,0,0);
 72:   if (tj->usehistory) {
 73:     TSHistoryUpdate(tj->tsh,stepnum,time);
 74:   }
 75:   if (tj->lag.caching) tj->lag.Udotcached.time = PETSC_MIN_REAL;
 76:   return(0);
 77: }

 79: /*@
 80:   TSTrajectoryGetNumSteps - Return the number of steps registered in the TSTrajectory via TSTrajectorySet().

 82:   Not collective.

 84:   Input Parameters:
 85: . tj - the trajectory object

 87:   Output Parameter:
 88: . steps - the number of steps

 90:   Level: developer

 92: .seealso: TSTrajectorySet()
 93: @*/
 94: PetscErrorCode TSTrajectoryGetNumSteps(TSTrajectory tj, PetscInt *steps)
 95: {

101:   TSHistoryGetNumSteps(tj->tsh,steps);
102:   return(0);
103: }

105: /*@
106:   TSTrajectoryGet - Updates the solution vector of a time stepper object by inquiring the TSTrajectory

108:   Collective on TS

110:   Input Parameters:
111: + tj      - the trajectory object
112: . ts      - the time stepper object
113: - stepnum - the step number

115:   Output Parameter:
116: . time    - the time associated with the step number

118:   Level: developer

120:   Notes: Usually one does not call this routine, it is called automatically during TSSolve()

122: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGetVecs(), TSGetSolution()
123: @*/
124: PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time)
125: {

129:   if (!tj) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS solver did not save trajectory");
134:   if (!tj->ops->get) SETERRQ1(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"TSTrajectory type %s",((PetscObject)tj)->type_name);
135:   if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
136:   if (stepnum < 0) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"Requesting negative step number");
137:   if (tj->monitor) {
138:     PetscViewerASCIIPrintf(tj->monitor,"TSTrajectoryGet: stepnum %D, stages %D\n",stepnum,(PetscInt)!tj->solution_only);
139:     PetscViewerFlush(tj->monitor);
140:   }
141:   PetscLogEventBegin(TSTrajectory_Get,tj,ts,0,0);
142:   (*tj->ops->get)(tj,ts,stepnum,time);
143:   PetscLogEventEnd(TSTrajectory_Get,tj,ts,0,0);
144:   return(0);
145: }

147: /*@
148:   TSTrajectoryGetVecs - Reconstructs the vector of state and its time derivative using information from the TSTrajectory and, possibly, from the TS

150:   Collective on TS

152:   Input Parameters:
153: + tj      - the trajectory object
154: . ts      - the time stepper object (optional)
155: - stepnum - the requested step number

157:   Input/Output Parameters:
158: . time - the time associated with the step number

160:   Output Parameters:
161: + U       - state vector (can be NULL)
162: - Udot    - time derivative of state vector (can be NULL)

164:   Level: developer

166:   Notes: If the step number is PETSC_DECIDE, the time argument is used to inquire the trajectory.
167:          If the requested time does not match any in the trajectory, Lagrangian interpolations are returned.

169: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGet()
170: @*/
171: PetscErrorCode TSTrajectoryGetVecs(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time,Vec U,Vec Udot)
172: {

176:   if (!tj) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS solver did not save trajectory");
183:   if (!U && !Udot) return(0);
184:   if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
185:   PetscLogEventBegin(TSTrajectory_GetVecs,tj,ts,0,0);
186:   if (tj->monitor) {
187:     PetscInt pU,pUdot;
188:     pU    = U ? 1 : 0;
189:     pUdot = Udot ? 1 : 0;
190:     PetscViewerASCIIPrintf(tj->monitor,"Requested by GetVecs %D %D: stepnum %D, time %g\n",pU,pUdot,stepnum,(double)*time);
191:     PetscViewerFlush(tj->monitor);
192:   }
193:   if (U && tj->lag.caching) {
194:     PetscObjectId    id;
195:     PetscObjectState state;

197:     PetscObjectStateGet((PetscObject)U,&state);
198:     PetscObjectGetId((PetscObject)U,&id);
199:     if (stepnum == PETSC_DECIDE) {
200:       if (id == tj->lag.Ucached.id && *time == tj->lag.Ucached.time && state == tj->lag.Ucached.state) U = NULL;
201:     } else {
202:       if (id == tj->lag.Ucached.id && stepnum == tj->lag.Ucached.step && state == tj->lag.Ucached.state) U = NULL;
203:     }
204:     if (tj->monitor && !U) {
205:       PetscViewerASCIIPushTab(tj->monitor);
206:       PetscViewerASCIIPrintf(tj->monitor,"State vector cached\n");
207:       PetscViewerASCIIPopTab(tj->monitor);
208:       PetscViewerFlush(tj->monitor);
209:     }
210:   }
211:   if (Udot && tj->lag.caching) {
212:     PetscObjectId    id;
213:     PetscObjectState state;

215:     PetscObjectStateGet((PetscObject)Udot,&state);
216:     PetscObjectGetId((PetscObject)Udot,&id);
217:     if (stepnum == PETSC_DECIDE) {
218:       if (id == tj->lag.Udotcached.id && *time == tj->lag.Udotcached.time && state == tj->lag.Udotcached.state) Udot = NULL;
219:     } else {
220:       if (id == tj->lag.Udotcached.id && stepnum == tj->lag.Udotcached.step && state == tj->lag.Udotcached.state) Udot = NULL;
221:     }
222:     if (tj->monitor && !Udot) {
223:       PetscViewerASCIIPushTab(tj->monitor);
224:       PetscViewerASCIIPrintf(tj->monitor,"Derivative vector cached\n");
225:       PetscViewerASCIIPopTab(tj->monitor);
226:       PetscViewerFlush(tj->monitor);
227:     }
228:   }
229:   if (!U && !Udot) {
230:     PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
231:     return(0);
232:   }

234:   if (stepnum == PETSC_DECIDE || Udot) { /* reverse search for requested time in TSHistory */
235:     if (tj->monitor) {
236:       PetscViewerASCIIPushTab(tj->monitor);
237:     }
238:     /* cached states will be updated in the function */
239:     TSTrajectoryReconstruct_Private(tj,ts,*time,U,Udot);
240:     if (tj->monitor) {
241:       PetscViewerASCIIPopTab(tj->monitor);
242:       PetscViewerFlush(tj->monitor);
243:     }
244:   } else if (U) { /* we were asked to load from stepnum, use TSTrajectoryGet */
245:     TS  fakets = ts;
246:     Vec U2;

248:     /* use a fake TS if ts is missing */
249:     if (!ts) {
250:       PetscObjectQuery((PetscObject)tj,"__fake_ts",(PetscObject*)&fakets);
251:       if (!fakets) {
252:         TSCreate(PetscObjectComm((PetscObject)tj),&fakets);
253:         PetscObjectCompose((PetscObject)tj,"__fake_ts",(PetscObject)fakets);
254:         PetscObjectDereference((PetscObject)fakets);
255:         VecDuplicate(U,&U2);
256:         TSSetSolution(fakets,U2);
257:         PetscObjectDereference((PetscObject)U2);
258:       }
259:     }
260:     TSTrajectoryGet(tj,fakets,stepnum,time);
261:     TSGetSolution(fakets,&U2);
262:     VecCopy(U2,U);
263:     PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
264:     PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
265:     tj->lag.Ucached.time = *time;
266:     tj->lag.Ucached.step = stepnum;
267:   }
268:   PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
269:   return(0);
270: }

272: /*@C
273:    TSTrajectoryViewFromOptions - View from Options

275:    Collective on TSTrajectory

277:    Input Parameters:
278: +  A - the TSTrajectory context
279: .  obj - Optional object
280: -  name - command line option

282:    Level: intermediate
283: .seealso:  TSTrajectory, TSTrajectoryView, PetscObjectViewFromOptions(), TSTrajectoryCreate()
284: @*/
285: PetscErrorCode  TSTrajectoryViewFromOptions(TSTrajectory A,PetscObject obj,const char name[])
286: {

291:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
292:   return(0);
293: }

295: /*@C
296:     TSTrajectoryView - Prints information about the trajectory object

298:     Collective on TSTrajectory

300:     Input Parameters:
301: +   tj - the TSTrajectory context obtained from TSTrajectoryCreate()
302: -   viewer - visualization context

304:     Options Database Key:
305: .   -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()

307:     Notes:
308:     The available visualization contexts include
309: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
310: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
311:          output where only the first processor opens
312:          the file.  All other processors send their
313:          data to the first processor to print.

315:     The user can open an alternative visualization context with
316:     PetscViewerASCIIOpen() - output to a specified file.

318:     Level: developer

320: .seealso: PetscViewerASCIIOpen()
321: @*/
322: PetscErrorCode  TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
323: {
325:   PetscBool      iascii;

329:   if (!viewer) {
330:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
331:   }

335:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
336:   if (iascii) {
337:     PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
338:     PetscViewerASCIIPrintf(viewer,"  total number of recomputations for adjoint calculation = %D\n",tj->recomps);
339:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint reads = %D\n",tj->diskreads);
340:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint writes = %D\n",tj->diskwrites);
341:     if (tj->ops->view) {
342:       PetscViewerASCIIPushTab(viewer);
343:       (*tj->ops->view)(tj,viewer);
344:       PetscViewerASCIIPopTab(viewer);
345:     }
346:   }
347:   return(0);
348: }

350: /*@C
351:    TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory

353:    Collective on TSTrajectory

355:    Input Parameters:
356: +  tr - the trajectory context
357: -  names - the names of the components, final string must be NULL

359:    Level: intermediate

361:    Note: Fortran interface is not possible because of the string array argument

363: .seealso: TSTrajectory, TSGetTrajectory()
364: @*/
365: PetscErrorCode  TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
366: {
367:   PetscErrorCode    ierr;

372:   PetscStrArrayDestroy(&ctx->names);
373:   PetscStrArrayallocpy(names,&ctx->names);
374:   return(0);
375: }

377: /*@C
378:    TSTrajectorySetTransform - Solution vector will be transformed by provided function before being saved to disk

380:    Collective on TSLGCtx

382:    Input Parameters:
383: +  tj - the TSTrajectory context
384: .  transform - the transform function
385: .  destroy - function to destroy the optional context
386: -  ctx - optional context used by transform function

388:    Level: intermediate

390: .seealso:  TSTrajectorySetVariableNames(), TSTrajectory, TSMonitorLGSetTransform()
391: @*/
392: PetscErrorCode  TSTrajectorySetTransform(TSTrajectory tj,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
393: {
396:   tj->transform        = transform;
397:   tj->transformdestroy = destroy;
398:   tj->transformctx     = tctx;
399:   return(0);
400: }

402: /*@
403:   TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE

405:   Collective

407:   Input Parameter:
408: . comm - the communicator

410:   Output Parameter:
411: . tj   - the trajectory object

413:   Level: developer

415:   Notes:
416:     Usually one does not call this routine, it is called automatically when one calls TSSetSaveTrajectory().

418: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
419: @*/
420: PetscErrorCode  TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
421: {
422:   TSTrajectory   t;

427:   *tj = NULL;
428:   TSInitializePackage();

430:   PetscHeaderCreate(t,TSTRAJECTORY_CLASSID,"TSTrajectory","Time stepping","TS",comm,TSTrajectoryDestroy,TSTrajectoryView);
431:   t->setupcalled = PETSC_FALSE;
432:   TSHistoryCreate(comm,&t->tsh);

434:   t->lag.order            = 1;
435:   t->lag.L                = NULL;
436:   t->lag.T                = NULL;
437:   t->lag.W                = NULL;
438:   t->lag.WW               = NULL;
439:   t->lag.TW               = NULL;
440:   t->lag.TT               = NULL;
441:   t->lag.caching          = PETSC_TRUE;
442:   t->lag.Ucached.id       = 0;
443:   t->lag.Ucached.state    = -1;
444:   t->lag.Ucached.time     = PETSC_MIN_REAL;
445:   t->lag.Ucached.step     = PETSC_MAX_INT;
446:   t->lag.Udotcached.id    = 0;
447:   t->lag.Udotcached.state = -1;
448:   t->lag.Udotcached.time  = PETSC_MIN_REAL;
449:   t->lag.Udotcached.step  = PETSC_MAX_INT;
450:   t->adjoint_solve_mode   = PETSC_TRUE;
451:   t->solution_only        = PETSC_FALSE;
452:   t->keepfiles            = PETSC_FALSE;
453:   t->usehistory           = PETSC_TRUE;
454:   *tj  = t;
455:   TSTrajectorySetFiletemplate(t,"TS-%06D.bin");
456:   return(0);
457: }

459: /*@C
460:   TSTrajectorySetType - Sets the storage method to be used as in a trajectory

462:   Collective on TS

464:   Input Parameters:
465: + tj   - the TSTrajectory context
466: . ts   - the TS context
467: - type - a known method

469:   Options Database Command:
470: . -ts_trajectory_type <type> - Sets the method; use -help for a list of available methods (for instance, basic)

472:    Level: developer

474: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy(), TSTrajectoryGetType()

476: @*/
477: PetscErrorCode  TSTrajectorySetType(TSTrajectory tj,TS ts,TSTrajectoryType type)
478: {
479:   PetscErrorCode (*r)(TSTrajectory,TS);
480:   PetscBool      match;

485:   PetscObjectTypeCompare((PetscObject)tj,type,&match);
486:   if (match) return(0);

488:   PetscFunctionListFind(TSTrajectoryList,type,&r);
489:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSTrajectory type: %s",type);
490:   if (tj->ops->destroy) {
491:     (*(tj)->ops->destroy)(tj);

493:     tj->ops->destroy = NULL;
494:   }
495:   PetscMemzero(tj->ops,sizeof(*tj->ops));

497:   PetscObjectChangeTypeName((PetscObject)tj,type);
498:   (*r)(tj,ts);
499:   return(0);
500: }

502: /*@C
503:   TSTrajectoryGetType - Gets the trajectory type

505:   Collective on TS

507:   Input Parameters:
508: + tj   - the TSTrajectory context
509: - ts   - the TS context

511:   Output Parameters:
512: . type - a known method

514:   Level: developer

516: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy(), TSTrajectorySetType()

518: @*/
519: PetscErrorCode TSTrajectoryGetType(TSTrajectory tj,TS ts,TSTrajectoryType *type)
520: {
523:   if (type) *type = ((PetscObject)tj)->type_name;
524:   return(0);
525: }

527: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
528: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
529: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
530: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);

532: /*@C
533:   TSTrajectoryRegisterAll - Registers all of the trajectory storage schecmes in the TS package.

535:   Not Collective

537:   Level: developer

539: .seealso: TSTrajectoryRegister()
540: @*/
541: PetscErrorCode  TSTrajectoryRegisterAll(void)
542: {

546:   if (TSTrajectoryRegisterAllCalled) return(0);
547:   TSTrajectoryRegisterAllCalled = PETSC_TRUE;

549:   TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
550:   TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
551:   TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
552:   TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
553:   return(0);
554: }

556: /*@
557:    TSTrajectoryReset - Resets a trajectory context

559:    Collective on TSTrajectory

561:    Input Parameter:
562: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

564:    Level: developer

566: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
567: @*/
568: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
569: {

573:   if (!tj) return(0);
575:   if (tj->ops->reset) {
576:     (*tj->ops->reset)(tj);
577:   }
578:   PetscFree(tj->dirfiletemplate);
579:   TSHistoryDestroy(&tj->tsh);
580:   TSHistoryCreate(PetscObjectComm((PetscObject)tj),&tj->tsh);
581:   tj->setupcalled = PETSC_FALSE;
582:   return(0);
583: }

585: /*@
586:    TSTrajectoryDestroy - Destroys a trajectory context

588:    Collective on TSTrajectory

590:    Input Parameter:
591: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

593:    Level: developer

595: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
596: @*/
597: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
598: {

602:   if (!*tj) return(0);
604:   if (--((PetscObject)(*tj))->refct > 0) {*tj = 0; return(0);}

606:   TSTrajectoryReset(*tj);
607:   TSHistoryDestroy(&(*tj)->tsh);
608:   VecDestroyVecs((*tj)->lag.order+1,&(*tj)->lag.W);
609:   PetscFree5((*tj)->lag.L,(*tj)->lag.T,(*tj)->lag.WW,(*tj)->lag.TT,(*tj)->lag.TW);
610:   VecDestroy(&(*tj)->U);
611:   VecDestroy(&(*tj)->Udot);

613:   if ((*tj)->transformdestroy) {(*(*tj)->transformdestroy)((*tj)->transformctx);}
614:   if ((*tj)->ops->destroy) {(*(*tj)->ops->destroy)((*tj));}
615:   if (!((*tj)->keepfiles)) {
616:     PetscMPIInt rank;
617:     MPI_Comm    comm;

619:     PetscObjectGetComm((PetscObject)(*tj),&comm);
620:     MPI_Comm_rank(comm,&rank);
621:     if (!rank && (*tj)->dirname) { /* we own the directory, so we run PetscRMTree on it */
622:       PetscRMTree((*tj)->dirname);
623:     }
624:   }
625:   PetscStrArrayDestroy(&(*tj)->names);
626:   PetscFree((*tj)->dirname);
627:   PetscFree((*tj)->filetemplate);
628:   PetscHeaderDestroy(tj);
629:   return(0);
630: }

632: /*
633:   TSTrajectorySetTypeFromOptions_Private - Sets the type of ts from user options.

635:   Collective on TSTrajectory

637:   Input Parameter:
638: + tj - the TSTrajectory context
639: - ts - the TS context

641:   Options Database Keys:
642: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION

644:   Level: developer

646: .seealso: TSTrajectorySetFromOptions(), TSTrajectorySetType()
647: */
648: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,TSTrajectory tj,TS ts)
649: {
650:   PetscBool      opt;
651:   const char     *defaultType;
652:   char           typeName[256];

656:   if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
657:   else defaultType = TSTRAJECTORYBASIC;

659:   TSTrajectoryRegisterAll();
660:   PetscOptionsFList("-ts_trajectory_type","TSTrajectory method","TSTrajectorySetType",TSTrajectoryList,defaultType,typeName,256,&opt);
661:   if (opt) {
662:     TSTrajectorySetType(tj,ts,typeName);
663:   } else {
664:     TSTrajectorySetType(tj,ts,defaultType);
665:   }
666:   return(0);
667: }

669: /*@
670:    TSTrajectorySetUseHistory - Use TSHistory in TSTrajectory

672:    Collective on TSTrajectory

674:    Input Arguments:
675: +  tj - the TSTrajectory context
676: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

678:    Options Database Keys:
679: .  -ts_trajectory_use_history - have it use TSHistory

681:    Level: advanced

683: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
684: @*/
685: PetscErrorCode TSTrajectorySetUseHistory(TSTrajectory tj,PetscBool flg)
686: {
690:   tj->usehistory = flg;
691:   return(0);
692: }

694: /*@
695:    TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller

697:    Collective on TSTrajectory

699:    Input Arguments:
700: +  tj - the TSTrajectory context
701: -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable

703:    Options Database Keys:
704: .  -ts_trajectory_monitor - print TSTrajectory information

706:    Level: developer

708: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
709: @*/
710: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
711: {
715:   if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
716:   else tj->monitor = NULL;
717:   return(0);
718: }

720: /*@
721:    TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory

723:    Collective on TSTrajectory

725:    Input Arguments:
726: +  tj - the TSTrajectory context
727: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

729:    Options Database Keys:
730: .  -ts_trajectory_keep_files - have it keep the files

732:    Notes:
733:     By default the TSTrajectory used for adjoint computations, TSTRAJECTORYBASIC, removes the files it generates at the end of the run. This causes the files to be kept.

735:    Level: advanced

737: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
738: @*/
739: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
740: {
744:   tj->keepfiles = flg;
745:   return(0);
746: }

748: /*@C
749:    TSTrajectorySetDirname - Specify the name of the directory where disk checkpoints are stored.

751:    Collective on TSTrajectory

753:    Input Arguments:
754: +  tj      - the TSTrajectory context
755: -  dirname - the directory name

757:    Options Database Keys:
758: .  -ts_trajectory_dirname - set the directory name

760:    Notes:
761:     The final location of the files is determined by dirname/filetemplate where filetemplate was provided by TSTrajectorySetFiletemplate()

763:    Level: developer

765: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
766: @*/
767: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
768: {
770:   PetscBool      flg;

774:   PetscStrcmp(tj->dirname,dirname,&flg);
775:   if (!flg && tj->dirfiletemplate) {
776:     SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after TSTrajectory has been setup");
777:   }
778:   PetscFree(tj->dirname);
779:   PetscStrallocpy(dirname,&tj->dirname);
780:   return(0);
781: }

783: /*@C
784:    TSTrajectorySetFiletemplate - Specify the name template for the files storing checkpoints.

786:    Collective on TSTrajectory

788:    Input Arguments:
789: +  tj      - the TSTrajectory context
790: -  filetemplate - the template

792:    Options Database Keys:
793: .  -ts_trajectory_file_template - set the file name template

795:    Notes:
796:     The name template should be of the form, for example filename-%06D.bin It should not begin with a leading /

798:    The final location of the files is determined by dirname/filetemplate where dirname was provided by TSTrajectorySetDirname(). The %06D is replaced by the
799:    timestep counter

801:    Level: developer

803: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
804: @*/
805: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
806: {
808:   const char     *ptr,*ptr2;

812:   if (tj->dirfiletemplate) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set filetemplate after TSTrajectory has been setup");

814:   if (!filetemplate[0]) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
815:   /* Do some cursory validation of the input. */
816:   PetscStrstr(filetemplate,"%",(char**)&ptr);
817:   if (!ptr) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
818:   for (ptr++; ptr && *ptr; ptr++) {
819:     PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
820:     if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"Invalid file template argument to -ts_trajectory_file_template, should look like filename-%%06D.bin");
821:     if (ptr2) break;
822:   }
823:   PetscFree(tj->filetemplate);
824:   PetscStrallocpy(filetemplate,&tj->filetemplate);
825:   return(0);
826: }

828: /*@
829:    TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.

831:    Collective on TSTrajectory

833:    Input Parameter:
834: +  tj - the TSTrajectory context obtained from TSTrajectoryCreate()
835: -  ts - the TS context

837:    Options Database Keys:
838: +  -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
839: .  -ts_trajectory_keep_files <true,false> - keep the files generated by the code after the program ends. This is true by default for TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
840: -  -ts_trajectory_monitor - print TSTrajectory information

842:    Level: developer

844:    Notes:
845:     This is not normally called directly by users

847: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
848: @*/
849: PetscErrorCode  TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
850: {
851:   PetscBool      set,flg;
852:   char           dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];

858:   PetscObjectOptionsBegin((PetscObject)tj);
859:   TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
860:   PetscOptionsBool("-ts_trajectory_use_history","Turn on/off usage of TSHistory",NULL,tj->usehistory,&tj->usehistory,NULL);
861:   PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
862:   if (set) {TSTrajectorySetMonitor(tj,flg);}
863:   PetscOptionsInt("-ts_trajectory_reconstruction_order","Interpolation order for reconstruction",NULL,tj->lag.order,&tj->lag.order,NULL);
864:   PetscOptionsBool("-ts_trajectory_reconstruction_caching","Turn on/off caching of TSTrajectoryGetVecs input",NULL,tj->lag.caching,&tj->lag.caching,NULL);
865:   PetscOptionsBool("-ts_trajectory_adjointmode","Instruct the trajectory that will be used in a TSAdjointSolve()",NULL,tj->adjoint_solve_mode,&tj->adjoint_solve_mode,NULL);
866:   PetscOptionsBool("-ts_trajectory_solution_only","Checkpoint solution only","TSTrajectorySetSolutionOnly",tj->solution_only,&tj->solution_only,NULL);
867:   PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
868:   if (set) {TSTrajectorySetKeepFiles(tj,flg);}

870:   PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",0,dirname,PETSC_MAX_PATH_LEN-14,&set);
871:   if (set) {
872:     TSTrajectorySetDirname(tj,dirname);
873:   }

875:   PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",0,filetemplate,PETSC_MAX_PATH_LEN,&set);
876:   if (set) {
877:     TSTrajectorySetFiletemplate(tj,filetemplate);
878:   }

880:   /* Handle specific TSTrajectory options */
881:   if (tj->ops->setfromoptions) {
882:     (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
883:   }
884:   PetscOptionsEnd();
885:   return(0);
886: }

888: /*@
889:    TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
890:    of a TS trajectory.

892:    Collective on TS

894:    Input Parameter:
895: +  ts - the TS context obtained from TSCreate()
896: -  tj - the TS trajectory context

898:    Level: developer

900: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
901: @*/
902: PetscErrorCode  TSTrajectorySetUp(TSTrajectory tj,TS ts)
903: {
905:   size_t         s1,s2;

908:   if (!tj) return(0);
911:   if (tj->setupcalled) return(0);

913:   if (!((PetscObject)tj)->type_name) {
914:     TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
915:   }
916:   if (tj->ops->setup) {
917:     (*tj->ops->setup)(tj,ts);
918:   }

920:   tj->setupcalled = PETSC_TRUE;

922:   /* Set the counters to zero */
923:   tj->recomps    = 0;
924:   tj->diskreads  = 0;
925:   tj->diskwrites = 0;
926:   PetscStrlen(tj->dirname,&s1);
927:   PetscStrlen(tj->filetemplate,&s2);
928:   PetscFree(tj->dirfiletemplate);
929:   PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
930:   PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
931:   return(0);
932: }

934: /*@
935:    TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage also.

937:    Collective on TSTrajectory

939:    Input Parameter:
940: +  tj  - the TS trajectory context
941: -  flg - the boolean flag

943:    Level: developer

945: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
946: @*/
947: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
948: {
952:   tj->solution_only = solution_only;
953:   return(0);
954: }

956: /*@
957:    TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.

959:    Logically collective on TSTrajectory

961:    Input Parameter:
962: .  tj  - the TS trajectory context

964:    Output Parameter:
965: -  flg - the boolean flag

967:    Level: developer

969: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetSolutionOnly()
970: @*/
971: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj,PetscBool *solution_only)
972: {
976:   *solution_only = tj->solution_only;
977:   return(0);
978: }

980: /*@
981:    TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.

983:    Collective on TSTrajectory

985:    Input Parameter:
986: +  tj   - the TS trajectory context
987: .  ts   - the TS solver context
988: -  time - the requested time

990:    Output Parameter:
991: +  U    - state vector at given time (can be interpolated)
992: -  Udot - time-derivative vector at given time (can be interpolated)

994:    Level: developer

996:    Notes: The vectors are interpolated if time does not match any time step stored in the TSTrajectory(). Pass NULL to not request a vector.
997:           This function differs from TSTrajectoryGetVecs since the vectors obtained cannot be modified, and they need to be returned by
998:           calling TSTrajectoryRestoreUpdatedHistoryVecs().

1000: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
1001: @*/
1002: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
1003: {

1012:   if (U && !tj->U) {
1013:     DM dm;

1015:     TSGetDM(ts,&dm);
1016:     DMCreateGlobalVector(dm,&tj->U);
1017:   }
1018:   if (Udot && !tj->Udot) {
1019:     DM dm;

1021:     TSGetDM(ts,&dm);
1022:     DMCreateGlobalVector(dm,&tj->Udot);
1023:   }
1024:   TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
1025:   if (U) {
1026:     VecLockReadPush(tj->U);
1027:     *U   = tj->U;
1028:   }
1029:   if (Udot) {
1030:     VecLockReadPush(tj->Udot);
1031:     *Udot = tj->Udot;
1032:   }
1033:   return(0);
1034: }

1036: /*@
1037:    TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with TSTrajectoryGetUpdatedHistoryVecs().

1039:    Collective on TSTrajectory

1041:    Input Parameter:
1042: +  tj   - the TS trajectory context
1043: .  U    - state vector at given time (can be interpolated)
1044: -  Udot - time-derivative vector at given time (can be interpolated)

1046:    Level: developer

1048: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
1049: @*/
1050: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1051: {

1058:   if (U && *U != tj->U) SETERRQ(PetscObjectComm((PetscObject)*U),PETSC_ERR_USER,"U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1059:   if (Udot && *Udot != tj->Udot) SETERRQ(PetscObjectComm((PetscObject)*Udot),PETSC_ERR_USER,"Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1060:   if (U) {
1061:     VecLockReadPop(tj->U);
1062:     *U   = NULL;
1063:   }
1064:   if (Udot) {
1065:     VecLockReadPop(tj->Udot);
1066:     *Udot = NULL;
1067:   }
1068:   return(0);
1069: }