Actual source code: traj.c
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, TSTrajectory_SetUp;
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 Parameter:
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 = NULL; 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 == 0 && (*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 Parameters:
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 Parameters:
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 Parameters:
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 Parameters:
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 Parameters:
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 Parameters:
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",NULL,dirname,sizeof(dirname)-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",NULL,filetemplate,sizeof(filetemplate),&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 Parameters:
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: PetscLogEventBegin(TSTrajectory_SetUp,tj,ts,0,0);
914: if (!((PetscObject)tj)->type_name) {
915: TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
916: }
917: if (tj->ops->setup) {
918: (*tj->ops->setup)(tj,ts);
919: }
921: tj->setupcalled = PETSC_TRUE;
923: /* Set the counters to zero */
924: tj->recomps = 0;
925: tj->diskreads = 0;
926: tj->diskwrites = 0;
927: PetscStrlen(tj->dirname,&s1);
928: PetscStrlen(tj->filetemplate,&s2);
929: PetscFree(tj->dirfiletemplate);
930: PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
931: PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
932: PetscLogEventEnd(TSTrajectory_SetUp,tj,ts,0,0);
933: return(0);
934: }
936: /*@
937: TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage also.
939: Collective on TSTrajectory
941: Input Parameters:
942: + tj - the TS trajectory context
943: - flg - the boolean flag
945: Level: developer
947: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
948: @*/
949: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
950: {
954: tj->solution_only = solution_only;
955: return(0);
956: }
958: /*@
959: TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.
961: Logically collective on TSTrajectory
963: Input Parameter:
964: . tj - the TS trajectory context
966: Output Parameter:
967: . flg - the boolean flag
969: Level: developer
971: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetSolutionOnly()
972: @*/
973: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj,PetscBool *solution_only)
974: {
978: *solution_only = tj->solution_only;
979: return(0);
980: }
982: /*@
983: TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.
985: Collective on TSTrajectory
987: Input Parameters:
988: + tj - the TS trajectory context
989: . ts - the TS solver context
990: - time - the requested time
992: Output Parameters:
993: + U - state vector at given time (can be interpolated)
994: - Udot - time-derivative vector at given time (can be interpolated)
996: Level: developer
998: Notes: The vectors are interpolated if time does not match any time step stored in the TSTrajectory(). Pass NULL to not request a vector.
999: This function differs from TSTrajectoryGetVecs since the vectors obtained cannot be modified, and they need to be returned by
1000: calling TSTrajectoryRestoreUpdatedHistoryVecs().
1002: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
1003: @*/
1004: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
1005: {
1014: if (U && !tj->U) {
1015: DM dm;
1017: TSGetDM(ts,&dm);
1018: DMCreateGlobalVector(dm,&tj->U);
1019: }
1020: if (Udot && !tj->Udot) {
1021: DM dm;
1023: TSGetDM(ts,&dm);
1024: DMCreateGlobalVector(dm,&tj->Udot);
1025: }
1026: TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
1027: if (U) {
1028: VecLockReadPush(tj->U);
1029: *U = tj->U;
1030: }
1031: if (Udot) {
1032: VecLockReadPush(tj->Udot);
1033: *Udot = tj->Udot;
1034: }
1035: return(0);
1036: }
1038: /*@
1039: TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with TSTrajectoryGetUpdatedHistoryVecs().
1041: Collective on TSTrajectory
1043: Input Parameters:
1044: + tj - the TS trajectory context
1045: . U - state vector at given time (can be interpolated)
1046: - Udot - time-derivative vector at given time (can be interpolated)
1048: Level: developer
1050: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
1051: @*/
1052: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1053: {
1060: if (U && *U != tj->U) SETERRQ(PetscObjectComm((PetscObject)*U),PETSC_ERR_USER,"U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1061: if (Udot && *Udot != tj->Udot) SETERRQ(PetscObjectComm((PetscObject)*Udot),PETSC_ERR_USER,"Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1062: if (U) {
1063: VecLockReadPop(tj->U);
1064: *U = NULL;
1065: }
1066: if (Udot) {
1067: VecLockReadPop(tj->Udot);
1068: *Udot = NULL;
1069: }
1070: return(0);
1071: }