Actual source code: traj.c

petsc-3.11.3 2019-06-26
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: .keywords: TS, trajectory, timestep, register

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

 33:   PetscFunctionListAdd(&TSTrajectoryList,sname,function);
 34:   return(0);
 35: }

 37: /*@
 38:   TSTrajectorySet - Sets a vector of state in the trajectory object

 40:   Collective on TSTrajectory

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

 49:   Level: developer

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

 53: .keywords: TS, trajectory, create

 55: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectoryGet(), TSTrajectoryGetVecs()
 56: @*/
 57: PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
 58: {

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

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

 84:   Not collective.

 86:   Input Parameters:
 87: . tj - the trajectory object

 89:   Output Parameter:
 90: . steps - the number of steps

 92:   Level: developer

 94: .keywords: TS, trajectory, create

 96: .seealso: TSTrajectorySet()
 97: @*/
 98: PetscErrorCode TSTrajectoryGetNumSteps(TSTrajectory tj, PetscInt *steps)
 99: {

105:   TSHistoryGetNumSteps(tj->tsh,steps);
106:   return(0);
107: }

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

112:   Collective on TS

114:   Input Parameters:
115: + tj      - the trajectory object
116: . ts      - the time stepper object
117: - stepnum - the step number

119:   Output Parameter:
120: . time    - the time associated with the step number

122:   Level: developer

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

126: .keywords: TS, trajectory, create

128: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGetVecs(), TSGetSolution()
129: @*/
130: PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time)
131: {

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

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

156:   Collective on TS

158:   Input Parameters:
159: + tj      - the trajectory object
160: . ts      - the time stepper object (optional)
161: - stepnum - the requested step number

163:   Input/Output Parameters:
164: . time - the time associated with the step number

166:   Output Parameters:
167: + U       - state vector (can be NULL)
168: - Udot    - time derivative of state vector (can be NULL)

170:   Level: developer

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

175: .keywords: TS, trajectory, create

177: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGet()
178: @*/
179: PetscErrorCode TSTrajectoryGetVecs(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time,Vec U,Vec Udot)
180: {

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

205:     PetscObjectStateGet((PetscObject)U,&state);
206:     PetscObjectGetId((PetscObject)U,&id);
207:     if (stepnum == PETSC_DECIDE) {
208:       if (id == tj->lag.Ucached.id && *time == tj->lag.Ucached.time && state == tj->lag.Ucached.state) U = NULL;
209:     } else {
210:       if (id == tj->lag.Ucached.id && stepnum == tj->lag.Ucached.step && state == tj->lag.Ucached.state) U = NULL;
211:     }
212:     if (tj->monitor && !U) {
213:       PetscViewerASCIIPushTab(tj->monitor);
214:       PetscViewerASCIIPrintf(tj->monitor,"State vector cached\n");
215:       PetscViewerASCIIPopTab(tj->monitor);
216:       PetscViewerFlush(tj->monitor);
217:     }
218:   }
219:   if (Udot && tj->lag.caching) {
220:     PetscObjectId    id;
221:     PetscObjectState state;

223:     PetscObjectStateGet((PetscObject)Udot,&state);
224:     PetscObjectGetId((PetscObject)Udot,&id);
225:     if (stepnum == PETSC_DECIDE) {
226:       if (id == tj->lag.Udotcached.id && *time == tj->lag.Udotcached.time && state == tj->lag.Udotcached.state) Udot = NULL;
227:     } else {
228:       if (id == tj->lag.Udotcached.id && stepnum == tj->lag.Udotcached.step && state == tj->lag.Udotcached.state) Udot = NULL;
229:     }
230:     if (tj->monitor && !Udot) {
231:       PetscViewerASCIIPushTab(tj->monitor);
232:       PetscViewerASCIIPrintf(tj->monitor,"Derivative vector cached\n");
233:       PetscViewerASCIIPopTab(tj->monitor);
234:       PetscViewerFlush(tj->monitor);
235:     }
236:   }
237:   if (!U && !Udot) {
238:     PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
239:     return(0);
240:   }

242:   if (stepnum == PETSC_DECIDE || Udot) { /* reverse search for requested time in TSHistory */
243:     if (tj->monitor) {
244:       PetscViewerASCIIPushTab(tj->monitor);
245:     }
246:     /* cached states will be updated in the function */
247:     TSTrajectoryReconstruct_Private(tj,ts,*time,U,Udot);
248:     if (tj->monitor) {
249:       PetscViewerASCIIPopTab(tj->monitor);
250:       PetscViewerFlush(tj->monitor);
251:     }
252:   } else if (U) { /* we were asked to load from stepnum, use TSTrajectoryGet */
253:     TS  fakets = ts;
254:     Vec U2;

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

280: /*@C
281:     TSTrajectoryView - Prints information about the trajectory object

283:     Collective on TSTrajectory

285:     Input Parameters:
286: +   tj - the TSTrajectory context obtained from TSTrajectoryCreate()
287: -   viewer - visualization context

289:     Options Database Key:
290: .   -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()

292:     Notes:
293:     The available visualization contexts include
294: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
295: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
296:          output where only the first processor opens
297:          the file.  All other processors send their
298:          data to the first processor to print.

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

303:     Level: developer

305: .keywords: TS, trajectory, timestep, view

307: .seealso: PetscViewerASCIIOpen()
308: @*/
309: PetscErrorCode  TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
310: {
312:   PetscBool      iascii;

316:   if (!viewer) {
317:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
318:   }

322:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
323:   if (iascii) {
324:     PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
325:     PetscViewerASCIIPrintf(viewer,"  total number of recomputations for adjoint calculation = %D\n",tj->recomps);
326:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint reads = %D\n",tj->diskreads);
327:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint writes = %D\n",tj->diskwrites);
328:     if (tj->ops->view) {
329:       PetscViewerASCIIPushTab(viewer);
330:       (*tj->ops->view)(tj,viewer);
331:       PetscViewerASCIIPopTab(viewer);
332:     }
333:   }
334:   return(0);
335: }

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

340:    Collective on TSTrajectory

342:    Input Parameters:
343: +  tr - the trajectory context
344: -  names - the names of the components, final string must be NULL

346:    Level: intermediate

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

350: .keywords: TS, TSTrajectory, vector, monitor, view

352: .seealso: TSTrajectory, TSGetTrajectory()
353: @*/
354: PetscErrorCode  TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
355: {
356:   PetscErrorCode    ierr;

361:   PetscStrArrayDestroy(&ctx->names);
362:   PetscStrArrayallocpy(names,&ctx->names);
363:   return(0);
364: }

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

369:    Collective on TSLGCtx

371:    Input Parameters:
372: +  tj - the TSTrajectory context
373: .  transform - the transform function
374: .  destroy - function to destroy the optional context
375: -  ctx - optional context used by transform function

377:    Level: intermediate

379: .keywords: TSTrajectory,  vector, monitor, view

381: .seealso:  TSTrajectorySetVariableNames(), TSTrajectory, TSMonitorLGSetTransform()
382: @*/
383: PetscErrorCode  TSTrajectorySetTransform(TSTrajectory tj,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
384: {
387:   tj->transform        = transform;
388:   tj->transformdestroy = destroy;
389:   tj->transformctx     = tctx;
390:   return(0);
391: }

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

396:   Collective on MPI_Comm

398:   Input Parameter:
399: . comm - the communicator

401:   Output Parameter:
402: . tj   - the trajectory object

404:   Level: developer

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

409: .keywords: TS, trajectory, create

411: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
412: @*/
413: PetscErrorCode  TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
414: {
415:   TSTrajectory   t;

420:   *tj = NULL;
421:   TSInitializePackage();

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

427:   t->lag.order            = 1;
428:   t->lag.L                = NULL;
429:   t->lag.T                = NULL;
430:   t->lag.W                = NULL;
431:   t->lag.WW               = NULL;
432:   t->lag.TW               = NULL;
433:   t->lag.TT               = NULL;
434:   t->lag.caching          = PETSC_TRUE;
435:   t->lag.Ucached.id       = 0;
436:   t->lag.Ucached.state    = -1;
437:   t->lag.Ucached.time     = PETSC_MIN_REAL;
438:   t->lag.Ucached.step     = PETSC_MAX_INT;
439:   t->lag.Udotcached.id    = 0;
440:   t->lag.Udotcached.state = -1;
441:   t->lag.Udotcached.time  = PETSC_MIN_REAL;
442:   t->lag.Udotcached.step  = PETSC_MAX_INT;
443:   t->adjoint_solve_mode   = PETSC_TRUE;
444:   t->solution_only        = PETSC_FALSE;
445:   t->keepfiles            = PETSC_FALSE;
446:   *tj  = t;
447:   TSTrajectorySetDirname(t,"SA-data");
448:   TSTrajectorySetFiletemplate(t,"SA-%06D.bin");
449:   return(0);
450: }

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

455:   Collective on TS

457:   Input Parameters:
458: + tj   - the TSTrajectory context
459: . ts   - the TS context
460: - type - a known method

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

465:    Level: developer

467: .keywords: TS, trajectory, timestep, set, type

469: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy()

471: @*/
472: PetscErrorCode  TSTrajectorySetType(TSTrajectory tj,TS ts,TSTrajectoryType type)
473: {
474:   PetscErrorCode (*r)(TSTrajectory,TS);
475:   PetscBool      match;

480:   PetscObjectTypeCompare((PetscObject)tj,type,&match);
481:   if (match) return(0);

483:   PetscFunctionListFind(TSTrajectoryList,type,&r);
484:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSTrajectory type: %s",type);
485:   if (tj->ops->destroy) {
486:     (*(tj)->ops->destroy)(tj);

488:     tj->ops->destroy = NULL;
489:   }
490:   PetscMemzero(tj->ops,sizeof(*tj->ops));

492:   PetscObjectChangeTypeName((PetscObject)tj,type);
493:   (*r)(tj,ts);
494:   return(0);
495: }

497: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
498: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
499: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
500: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);

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

505:   Not Collective

507:   Level: developer

509: .keywords: TS, trajectory, register, all

511: .seealso: TSTrajectoryRegister()
512: @*/
513: PetscErrorCode  TSTrajectoryRegisterAll(void)
514: {

518:   if (TSTrajectoryRegisterAllCalled) return(0);
519:   TSTrajectoryRegisterAllCalled = PETSC_TRUE;

521:   TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
522:   TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
523:   TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
524:   TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
525:   return(0);
526: }

528: /*@
529:    TSTrajectoryReset - Resets a trajectory context

531:    Collective on TSTrajectory

533:    Input Parameter:
534: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

536:    Level: developer

538: .keywords: TS, trajectory, timestep, reset

540: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
541: @*/
542: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
543: {

547:   if (!tj) return(0);
549:   if (tj->ops->reset) {
550:     (*tj->ops->reset)(tj);
551:   }
552:   PetscFree(tj->dirfiletemplate);
553:   TSHistoryDestroy(&tj->tsh);
554:   TSHistoryCreate(PetscObjectComm((PetscObject)tj),&tj->tsh);
555:   tj->setupcalled = PETSC_FALSE;
556:   return(0);
557: }

559: /*@
560:    TSTrajectoryDestroy - Destroys a trajectory context

562:    Collective on TSTrajectory

564:    Input Parameter:
565: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

567:    Level: developer

569: .keywords: TS, trajectory, timestep, destroy

571: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
572: @*/
573: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
574: {

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

582:   TSTrajectoryReset(*tj);
583:   TSHistoryDestroy(&(*tj)->tsh);
584:   VecDestroyVecs((*tj)->lag.order+1,&(*tj)->lag.W);
585:   PetscFree5((*tj)->lag.L,(*tj)->lag.T,(*tj)->lag.WW,(*tj)->lag.TT,(*tj)->lag.TW);
586:   VecDestroy(&(*tj)->U);
587:   VecDestroy(&(*tj)->Udot);

589:   if ((*tj)->transformdestroy) {(*(*tj)->transformdestroy)((*tj)->transformctx);}
590:   if ((*tj)->ops->destroy) {(*(*tj)->ops->destroy)((*tj));}
591:   if (!((*tj)->keepfiles)) {
592:     PetscMPIInt rank;
593:     MPI_Comm    comm;

595:     PetscObjectGetComm((PetscObject)(*tj),&comm);
596:     MPI_Comm_rank(comm,&rank);
597:     if (!rank) { /* we own the directory, so we run PetscRMTree on it */
598:       PetscRMTree((*tj)->dirname);
599:     }
600:   }
601:   PetscStrArrayDestroy(&(*tj)->names);
602:   PetscFree((*tj)->dirname);
603:   PetscFree((*tj)->filetemplate);
604:   PetscHeaderDestroy(tj);
605:   return(0);
606: }

608: /*
609:   TSTrajectorySetTypeFromOptions_Private - Sets the type of ts from user options.

611:   Collective on TSTrajectory

613:   Input Parameter:
614: + tj - the TSTrajectory context
615: - ts - the TS context

617:   Options Database Keys:
618: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION

620:   Level: developer

622: .keywords: TS, trajectory, set, options, type

624: .seealso: TSTrajectorySetFromOptions(), TSTrajectorySetType()
625: */
626: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,TSTrajectory tj,TS ts)
627: {
628:   PetscBool      opt;
629:   const char     *defaultType;
630:   char           typeName[256];

634:   if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
635:   else defaultType = TSTRAJECTORYBASIC;

637:   TSTrajectoryRegisterAll();
638:   PetscOptionsFList("-ts_trajectory_type","TSTrajectory method","TSTrajectorySetType",TSTrajectoryList,defaultType,typeName,256,&opt);
639:   if (opt) {
640:     TSTrajectorySetType(tj,ts,typeName);
641:   } else {
642:     TSTrajectorySetType(tj,ts,defaultType);
643:   }
644:   return(0);
645: }

647: /*@
648:    TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller

650:    Collective on TSTrajectory

652:    Input Arguments:
653: +  tj - the TSTrajectory context
654: -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable

656:    Options Database Keys:
657: .  -ts_trajectory_monitor - print TSTrajectory information

659:    Level: developer

661: .keywords: TS, trajectory, set, monitor

663: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
664: @*/
665: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
666: {
670:   if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
671:   else tj->monitor = NULL;
672:   return(0);
673: }

675: /*@
676:    TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory

678:    Collective on TSTrajectory

680:    Input Arguments:
681: +  tj - the TSTrajectory context
682: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

684:    Options Database Keys:
685: .  -ts_trajectory_keep_files - have it keep the files

687:    Notes:
688:     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.

690:    Level: advanced

692: .keywords: TS, trajectory, set, monitor

694: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
695: @*/
696: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
697: {
701:   tj->keepfiles = flg;
702:   return(0);
703: }

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

708:    Collective on TSTrajectory

710:    Input Arguments:
711: +  tj      - the TSTrajectory context
712: -  dirname - the directory name

714:    Options Database Keys:
715: .  -ts_trajectory_dirname - set the directory name

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

720:    Level: developer

722: .keywords: TS, trajectory, set

724: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
725: @*/
726: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
727: {
729:   PetscBool      flg;

733:   PetscStrcmp(tj->dirname,dirname,&flg);
734:   if (!flg && tj->dirfiletemplate) {
735:     SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after TSTrajectory has been setup");
736:   }
737:   PetscFree(tj->dirname);
738:   PetscStrallocpy(dirname,&tj->dirname);
739:   return(0);
740: }

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

745:    Collective on TSTrajectory

747:    Input Arguments:
748: +  tj      - the TSTrajectory context
749: -  filetemplate - the template

751:    Options Database Keys:
752: .  -ts_trajectory_file_template - set the file name template

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

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

760:    Level: developer

762: .keywords: TS, trajectory, set

764: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
765: @*/
766: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
767: {
769:   const char     *ptr,*ptr2;

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

775:   if (!filetemplate[0]) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
776:   /* Do some cursory validation of the input. */
777:   PetscStrstr(filetemplate,"%",(char**)&ptr);
778:   if (!ptr) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
779:   for (ptr++; ptr && *ptr; ptr++) {
780:     PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
781:     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");
782:     if (ptr2) break;
783:   }
784:   PetscFree(tj->filetemplate);
785:   PetscStrallocpy(filetemplate,&tj->filetemplate);
786:   return(0);
787: }

789: /*@
790:    TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.

792:    Collective on TSTrajectory

794:    Input Parameter:
795: +  tj - the TSTrajectory context obtained from TSTrajectoryCreate()
796: -  ts - the TS context

798:    Options Database Keys:
799: +  -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
800: .  -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
801: -  -ts_trajectory_monitor - print TSTrajectory information

803:    Level: developer

805:    Notes:
806:     This is not normally called directly by users

808: .keywords: TS, trajectory, timestep, set, options, database

810: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
811: @*/
812: PetscErrorCode  TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
813: {
814:   PetscBool      set,flg;
815:   char           dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];

821:   PetscObjectOptionsBegin((PetscObject)tj);
822:   TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
823:   PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
824:   if (set) {TSTrajectorySetMonitor(tj,flg);}
825:   PetscOptionsInt("-ts_trajectory_reconstruction_order","Interpolation order for reconstruction",NULL,tj->lag.order,&tj->lag.order,NULL);
826:   PetscOptionsBool("-ts_trajectory_reconstruction_caching","Turn on/off caching of TSTrajectoryGetVecs input",NULL,tj->lag.caching,&tj->lag.caching,NULL);
827:   PetscOptionsBool("-ts_trajectory_adjointmode","Instruct the trajectory that will be used in a TSAdjointSolve()",NULL,tj->adjoint_solve_mode,&tj->adjoint_solve_mode,NULL);
828:   PetscOptionsBool("-ts_trajectory_solution_only","Checkpoint solution only","TSTrajectorySetSolutionOnly",tj->solution_only,&tj->solution_only,NULL);
829:   PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
830:   if (set) {TSTrajectorySetKeepFiles(tj,flg);}

832:   PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",0,dirname,PETSC_MAX_PATH_LEN-14,&set);
833:   if (set) {
834:     TSTrajectorySetDirname(tj,dirname);
835:   }

837:   PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",0,filetemplate,PETSC_MAX_PATH_LEN,&set);
838:   if (set) {
839:     TSTrajectorySetFiletemplate(tj,filetemplate);
840:   }

842:   /* Handle specific TSTrajectory options */
843:   if (tj->ops->setfromoptions) {
844:     (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
845:   }
846:   PetscOptionsEnd();
847:   return(0);
848: }

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

854:    Collective on TS

856:    Input Parameter:
857: +  ts - the TS context obtained from TSCreate()
858: -  tj - the TS trajectory context

860:    Level: developer

862: .keywords: TS, trajectory, setup

864: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
865: @*/
866: PetscErrorCode  TSTrajectorySetUp(TSTrajectory tj,TS ts)
867: {
869:   size_t         s1,s2;

872:   if (!tj) return(0);
875:   if (tj->setupcalled) return(0);

877:   if (!((PetscObject)tj)->type_name) {
878:     TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
879:   }
880:   if (tj->ops->setup) {
881:     (*tj->ops->setup)(tj,ts);
882:   }

884:   tj->setupcalled = PETSC_TRUE;

886:   /* Set the counters to zero */
887:   tj->recomps    = 0;
888:   tj->diskreads  = 0;
889:   tj->diskwrites = 0;
890:   PetscStrlen(tj->dirname,&s1);
891:   PetscStrlen(tj->filetemplate,&s2);
892:   PetscFree(tj->dirfiletemplate);
893:   PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
894:   PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
895:   return(0);
896: }

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

901:    Collective on TSTrajectory

903:    Input Parameter:
904: +  tj  - the TS trajectory context
905: -  flg - the boolean flag

907:    Level: developer

909: .keywords: trajectory

911: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
912: @*/
913: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
914: {
918:   tj->solution_only = solution_only;
919:   return(0);
920: }

922: /*@
923:    TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.

925:    Logically collective on TSTrajectory

927:    Input Parameter:
928: .  tj  - the TS trajectory context

930:    Output Parameter:
931: -  flg - the boolean flag

933:    Level: developer

935: .keywords: trajectory

937: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetSolutionOnly()
938: @*/
939: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj,PetscBool *solution_only)
940: {
944:   *solution_only = tj->solution_only;
945:   return(0);
946: }

948: /*@
949:    TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.

951:    Collective on TSTrajectory

953:    Input Parameter:
954: +  tj   - the TS trajectory context
955: .  ts   - the TS solver context
956: -  time - the requested time

958:    Output Parameter:
959: +  U    - state vector at given time (can be interpolated)
960: -  Udot - time-derivative vector at given time (can be interpolated)

962:    Level: developer

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

968: .keywords: trajectory

970: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
971: @*/
972: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
973: {

982:   if (U && !tj->U) {
983:     DM dm;

985:     TSGetDM(ts,&dm);
986:     DMCreateGlobalVector(dm,&tj->U);
987:   }
988:   if (Udot && !tj->Udot) {
989:     DM dm;

991:     TSGetDM(ts,&dm);
992:     DMCreateGlobalVector(dm,&tj->Udot);
993:   }
994:   TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
995:   if (U) {
996:     VecLockReadPush(tj->U);
997:     *U   = tj->U;
998:   }
999:   if (Udot) {
1000:     VecLockReadPush(tj->Udot);
1001:     *Udot = tj->Udot;
1002:   }
1003:   return(0);
1004: }

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

1009:    Collective on TSTrajectory

1011:    Input Parameter:
1012: +  tj   - the TS trajectory context
1013: .  U    - state vector at given time (can be interpolated)
1014: -  Udot - time-derivative vector at given time (can be interpolated)

1016:    Level: developer

1018: .keywords: trajectory

1020: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
1021: @*/
1022: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1023: {

1030:   if (U && *U != tj->U) SETERRQ(PetscObjectComm((PetscObject)*U),PETSC_ERR_USER,"U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1031:   if (Udot && *Udot != tj->Udot) SETERRQ(PetscObjectComm((PetscObject)*Udot),PETSC_ERR_USER,"Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1032:   if (U) {
1033:     VecLockReadPop(tj->U);
1034:     *U   = NULL;
1035:   }
1036:   if (Udot) {
1037:     VecLockReadPop(tj->Udot);
1038:     *Udot = NULL;
1039:   }
1040:   return(0);
1041: }