Actual source code: traj.c

petsc-master 2019-09-18
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:     TSTrajectoryView - Prints information about the trajectory object

275:     Collective on TSTrajectory

277:     Input Parameters:
278: +   tj - the TSTrajectory context obtained from TSTrajectoryCreate()
279: -   viewer - visualization context

281:     Options Database Key:
282: .   -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()

284:     Notes:
285:     The available visualization contexts include
286: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
287: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
288:          output where only the first processor opens
289:          the file.  All other processors send their
290:          data to the first processor to print.

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

295:     Level: developer

297: .seealso: PetscViewerASCIIOpen()
298: @*/
299: PetscErrorCode  TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
300: {
302:   PetscBool      iascii;

306:   if (!viewer) {
307:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
308:   }

312:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
313:   if (iascii) {
314:     PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
315:     PetscViewerASCIIPrintf(viewer,"  total number of recomputations for adjoint calculation = %D\n",tj->recomps);
316:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint reads = %D\n",tj->diskreads);
317:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint writes = %D\n",tj->diskwrites);
318:     if (tj->ops->view) {
319:       PetscViewerASCIIPushTab(viewer);
320:       (*tj->ops->view)(tj,viewer);
321:       PetscViewerASCIIPopTab(viewer);
322:     }
323:   }
324:   return(0);
325: }

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

330:    Collective on TSTrajectory

332:    Input Parameters:
333: +  tr - the trajectory context
334: -  names - the names of the components, final string must be NULL

336:    Level: intermediate

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

340: .seealso: TSTrajectory, TSGetTrajectory()
341: @*/
342: PetscErrorCode  TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
343: {
344:   PetscErrorCode    ierr;

349:   PetscStrArrayDestroy(&ctx->names);
350:   PetscStrArrayallocpy(names,&ctx->names);
351:   return(0);
352: }

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

357:    Collective on TSLGCtx

359:    Input Parameters:
360: +  tj - the TSTrajectory context
361: .  transform - the transform function
362: .  destroy - function to destroy the optional context
363: -  ctx - optional context used by transform function

365:    Level: intermediate

367: .seealso:  TSTrajectorySetVariableNames(), TSTrajectory, TSMonitorLGSetTransform()
368: @*/
369: PetscErrorCode  TSTrajectorySetTransform(TSTrajectory tj,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
370: {
373:   tj->transform        = transform;
374:   tj->transformdestroy = destroy;
375:   tj->transformctx     = tctx;
376:   return(0);
377: }

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

382:   Collective

384:   Input Parameter:
385: . comm - the communicator

387:   Output Parameter:
388: . tj   - the trajectory object

390:   Level: developer

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

395: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
396: @*/
397: PetscErrorCode  TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
398: {
399:   TSTrajectory   t;

404:   *tj = NULL;
405:   TSInitializePackage();

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

411:   t->lag.order            = 1;
412:   t->lag.L                = NULL;
413:   t->lag.T                = NULL;
414:   t->lag.W                = NULL;
415:   t->lag.WW               = NULL;
416:   t->lag.TW               = NULL;
417:   t->lag.TT               = NULL;
418:   t->lag.caching          = PETSC_TRUE;
419:   t->lag.Ucached.id       = 0;
420:   t->lag.Ucached.state    = -1;
421:   t->lag.Ucached.time     = PETSC_MIN_REAL;
422:   t->lag.Ucached.step     = PETSC_MAX_INT;
423:   t->lag.Udotcached.id    = 0;
424:   t->lag.Udotcached.state = -1;
425:   t->lag.Udotcached.time  = PETSC_MIN_REAL;
426:   t->lag.Udotcached.step  = PETSC_MAX_INT;
427:   t->adjoint_solve_mode   = PETSC_TRUE;
428:   t->solution_only        = PETSC_FALSE;
429:   t->keepfiles            = PETSC_FALSE;
430:   t->usehistory           = PETSC_TRUE;
431:   *tj  = t;
432:   TSTrajectorySetFiletemplate(t,"TS-%06D.bin");
433:   return(0);
434: }

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

439:   Collective on TS

441:   Input Parameters:
442: + tj   - the TSTrajectory context
443: . ts   - the TS context
444: - type - a known method

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

449:    Level: developer

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

453: @*/
454: PetscErrorCode  TSTrajectorySetType(TSTrajectory tj,TS ts,TSTrajectoryType type)
455: {
456:   PetscErrorCode (*r)(TSTrajectory,TS);
457:   PetscBool      match;

462:   PetscObjectTypeCompare((PetscObject)tj,type,&match);
463:   if (match) return(0);

465:   PetscFunctionListFind(TSTrajectoryList,type,&r);
466:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSTrajectory type: %s",type);
467:   if (tj->ops->destroy) {
468:     (*(tj)->ops->destroy)(tj);

470:     tj->ops->destroy = NULL;
471:   }
472:   PetscMemzero(tj->ops,sizeof(*tj->ops));

474:   PetscObjectChangeTypeName((PetscObject)tj,type);
475:   (*r)(tj,ts);
476:   return(0);
477: }

479: /*@C
480:   TSTrajectoryGetType - Gets the trajectory type

482:   Collective on TS

484:   Input Parameters:
485: + tj   - the TSTrajectory context
486: - ts   - the TS context

488:   Output Parameters:
489: . type - a known method

491:   Level: developer

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

495: @*/
496: PetscErrorCode TSTrajectoryGetType(TSTrajectory tj,TS ts,TSTrajectoryType *type)
497: {
500:   if (type) *type = ((PetscObject)tj)->type_name;
501:   return(0);
502: }

504: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
505: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
506: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
507: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);

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

512:   Not Collective

514:   Level: developer

516: .seealso: TSTrajectoryRegister()
517: @*/
518: PetscErrorCode  TSTrajectoryRegisterAll(void)
519: {

523:   if (TSTrajectoryRegisterAllCalled) return(0);
524:   TSTrajectoryRegisterAllCalled = PETSC_TRUE;

526:   TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
527:   TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
528:   TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
529:   TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
530:   return(0);
531: }

533: /*@
534:    TSTrajectoryReset - Resets a trajectory context

536:    Collective on TSTrajectory

538:    Input Parameter:
539: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

541:    Level: developer

543: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
544: @*/
545: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
546: {

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

562: /*@
563:    TSTrajectoryDestroy - Destroys a trajectory context

565:    Collective on TSTrajectory

567:    Input Parameter:
568: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

570:    Level: developer

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

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

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

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

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

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

612:   Collective on TSTrajectory

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

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

621:   Level: developer

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

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

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

646: /*@
647:    TSTrajectorySetUseHistory - Use TSHistory in TSTrajectory

649:    Collective on TSTrajectory

651:    Input Arguments:
652: +  tj - the TSTrajectory context
653: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

655:    Options Database Keys:
656: .  -ts_trajectory_use_history - have it use TSHistory

658:    Level: advanced

660: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
661: @*/
662: PetscErrorCode TSTrajectorySetUseHistory(TSTrajectory tj,PetscBool flg)
663: {
667:   tj->usehistory = flg;
668:   return(0);
669: }

671: /*@
672:    TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller

674:    Collective on TSTrajectory

676:    Input Arguments:
677: +  tj - the TSTrajectory context
678: -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable

680:    Options Database Keys:
681: .  -ts_trajectory_monitor - print TSTrajectory information

683:    Level: developer

685: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
686: @*/
687: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
688: {
692:   if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
693:   else tj->monitor = NULL;
694:   return(0);
695: }

697: /*@
698:    TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory

700:    Collective on TSTrajectory

702:    Input Arguments:
703: +  tj - the TSTrajectory context
704: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

706:    Options Database Keys:
707: .  -ts_trajectory_keep_files - have it keep the files

709:    Notes:
710:     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.

712:    Level: advanced

714: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
715: @*/
716: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
717: {
721:   tj->keepfiles = flg;
722:   return(0);
723: }

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

728:    Collective on TSTrajectory

730:    Input Arguments:
731: +  tj      - the TSTrajectory context
732: -  dirname - the directory name

734:    Options Database Keys:
735: .  -ts_trajectory_dirname - set the directory name

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

740:    Level: developer

742: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
743: @*/
744: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
745: {
747:   PetscBool      flg;

751:   PetscStrcmp(tj->dirname,dirname,&flg);
752:   if (!flg && tj->dirfiletemplate) {
753:     SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after TSTrajectory has been setup");
754:   }
755:   PetscFree(tj->dirname);
756:   PetscStrallocpy(dirname,&tj->dirname);
757:   return(0);
758: }

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

763:    Collective on TSTrajectory

765:    Input Arguments:
766: +  tj      - the TSTrajectory context
767: -  filetemplate - the template

769:    Options Database Keys:
770: .  -ts_trajectory_file_template - set the file name template

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

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

778:    Level: developer

780: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
781: @*/
782: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
783: {
785:   const char     *ptr,*ptr2;

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

791:   if (!filetemplate[0]) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
792:   /* Do some cursory validation of the input. */
793:   PetscStrstr(filetemplate,"%",(char**)&ptr);
794:   if (!ptr) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
795:   for (ptr++; ptr && *ptr; ptr++) {
796:     PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
797:     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");
798:     if (ptr2) break;
799:   }
800:   PetscFree(tj->filetemplate);
801:   PetscStrallocpy(filetemplate,&tj->filetemplate);
802:   return(0);
803: }

805: /*@
806:    TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.

808:    Collective on TSTrajectory

810:    Input Parameter:
811: +  tj - the TSTrajectory context obtained from TSTrajectoryCreate()
812: -  ts - the TS context

814:    Options Database Keys:
815: +  -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
816: .  -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
817: -  -ts_trajectory_monitor - print TSTrajectory information

819:    Level: developer

821:    Notes:
822:     This is not normally called directly by users

824: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
825: @*/
826: PetscErrorCode  TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
827: {
828:   PetscBool      set,flg;
829:   char           dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];

835:   PetscObjectOptionsBegin((PetscObject)tj);
836:   TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
837:   PetscOptionsBool("-ts_trajectory_use_history","Turn on/off usage of TSHistory",NULL,tj->usehistory,&tj->usehistory,NULL);
838:   PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
839:   if (set) {TSTrajectorySetMonitor(tj,flg);}
840:   PetscOptionsInt("-ts_trajectory_reconstruction_order","Interpolation order for reconstruction",NULL,tj->lag.order,&tj->lag.order,NULL);
841:   PetscOptionsBool("-ts_trajectory_reconstruction_caching","Turn on/off caching of TSTrajectoryGetVecs input",NULL,tj->lag.caching,&tj->lag.caching,NULL);
842:   PetscOptionsBool("-ts_trajectory_adjointmode","Instruct the trajectory that will be used in a TSAdjointSolve()",NULL,tj->adjoint_solve_mode,&tj->adjoint_solve_mode,NULL);
843:   PetscOptionsBool("-ts_trajectory_solution_only","Checkpoint solution only","TSTrajectorySetSolutionOnly",tj->solution_only,&tj->solution_only,NULL);
844:   PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
845:   if (set) {TSTrajectorySetKeepFiles(tj,flg);}

847:   PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",0,dirname,PETSC_MAX_PATH_LEN-14,&set);
848:   if (set) {
849:     TSTrajectorySetDirname(tj,dirname);
850:   }

852:   PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",0,filetemplate,PETSC_MAX_PATH_LEN,&set);
853:   if (set) {
854:     TSTrajectorySetFiletemplate(tj,filetemplate);
855:   }

857:   /* Handle specific TSTrajectory options */
858:   if (tj->ops->setfromoptions) {
859:     (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
860:   }
861:   PetscOptionsEnd();
862:   return(0);
863: }

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

869:    Collective on TS

871:    Input Parameter:
872: +  ts - the TS context obtained from TSCreate()
873: -  tj - the TS trajectory context

875:    Level: developer

877: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
878: @*/
879: PetscErrorCode  TSTrajectorySetUp(TSTrajectory tj,TS ts)
880: {
882:   size_t         s1,s2;

885:   if (!tj) return(0);
888:   if (tj->setupcalled) return(0);

890:   if (!((PetscObject)tj)->type_name) {
891:     TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
892:   }
893:   if (tj->ops->setup) {
894:     (*tj->ops->setup)(tj,ts);
895:   }

897:   tj->setupcalled = PETSC_TRUE;

899:   /* Set the counters to zero */
900:   tj->recomps    = 0;
901:   tj->diskreads  = 0;
902:   tj->diskwrites = 0;
903:   PetscStrlen(tj->dirname,&s1);
904:   PetscStrlen(tj->filetemplate,&s2);
905:   PetscFree(tj->dirfiletemplate);
906:   PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
907:   PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
908:   return(0);
909: }

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

914:    Collective on TSTrajectory

916:    Input Parameter:
917: +  tj  - the TS trajectory context
918: -  flg - the boolean flag

920:    Level: developer

922: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
923: @*/
924: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
925: {
929:   tj->solution_only = solution_only;
930:   return(0);
931: }

933: /*@
934:    TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.

936:    Logically collective on TSTrajectory

938:    Input Parameter:
939: .  tj  - the TS trajectory context

941:    Output Parameter:
942: -  flg - the boolean flag

944:    Level: developer

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

957: /*@
958:    TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.

960:    Collective on TSTrajectory

962:    Input Parameter:
963: +  tj   - the TS trajectory context
964: .  ts   - the TS solver context
965: -  time - the requested time

967:    Output Parameter:
968: +  U    - state vector at given time (can be interpolated)
969: -  Udot - time-derivative vector at given time (can be interpolated)

971:    Level: developer

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

977: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
978: @*/
979: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
980: {

989:   if (U && !tj->U) {
990:     DM dm;

992:     TSGetDM(ts,&dm);
993:     DMCreateGlobalVector(dm,&tj->U);
994:   }
995:   if (Udot && !tj->Udot) {
996:     DM dm;

998:     TSGetDM(ts,&dm);
999:     DMCreateGlobalVector(dm,&tj->Udot);
1000:   }
1001:   TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
1002:   if (U) {
1003:     VecLockReadPush(tj->U);
1004:     *U   = tj->U;
1005:   }
1006:   if (Udot) {
1007:     VecLockReadPush(tj->Udot);
1008:     *Udot = tj->Udot;
1009:   }
1010:   return(0);
1011: }

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

1016:    Collective on TSTrajectory

1018:    Input Parameter:
1019: +  tj   - the TS trajectory context
1020: .  U    - state vector at given time (can be interpolated)
1021: -  Udot - time-derivative vector at given time (can be interpolated)

1023:    Level: developer

1025: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
1026: @*/
1027: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1028: {

1035:   if (U && *U != tj->U) SETERRQ(PetscObjectComm((PetscObject)*U),PETSC_ERR_USER,"U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1036:   if (Udot && *Udot != tj->Udot) SETERRQ(PetscObjectComm((PetscObject)*Udot),PETSC_ERR_USER,"Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1037:   if (U) {
1038:     VecLockReadPop(tj->U);
1039:     *U   = NULL;
1040:   }
1041:   if (Udot) {
1042:     VecLockReadPop(tj->Udot);
1043:     *Udot = NULL;
1044:   }
1045:   return(0);
1046: }