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: {
28: PetscFunctionListAdd(&TSTrajectoryList,sname,function);
29: return 0;
30: }
32: /*@
33: TSTrajectorySet - Sets a vector of state in the trajectory object
35: Collective on TSTrajectory
37: Input Parameters:
38: + tj - the trajectory object
39: . ts - the time stepper object (optional)
40: . stepnum - the step number
41: . time - the current time
42: - X - the current solution
44: Level: developer
46: Notes: Usually one does not call this routine, it is called automatically during TSSolve()
48: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectoryGet(), TSTrajectoryGetVecs()
49: @*/
50: PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
51: {
52: if (!tj) return 0;
60: if (tj->monitor) {
61: PetscViewerASCIIPrintf(tj->monitor,"TSTrajectorySet: stepnum %D, time %g (stages %D)\n",stepnum,(double)time,(PetscInt)!tj->solution_only);
62: }
63: PetscLogEventBegin(TSTrajectory_Set,tj,ts,0,0);
64: (*tj->ops->set)(tj,ts,stepnum,time,X);
65: PetscLogEventEnd(TSTrajectory_Set,tj,ts,0,0);
66: if (tj->usehistory) {
67: TSHistoryUpdate(tj->tsh,stepnum,time);
68: }
69: if (tj->lag.caching) tj->lag.Udotcached.time = PETSC_MIN_REAL;
70: return 0;
71: }
73: /*@
74: TSTrajectoryGetNumSteps - Return the number of steps registered in the TSTrajectory via TSTrajectorySet().
76: Not collective.
78: Input Parameters:
79: . tj - the trajectory object
81: Output Parameter:
82: . steps - the number of steps
84: Level: developer
86: .seealso: TSTrajectorySet()
87: @*/
88: PetscErrorCode TSTrajectoryGetNumSteps(TSTrajectory tj, PetscInt *steps)
89: {
92: TSHistoryGetNumSteps(tj->tsh,steps);
93: return 0;
94: }
96: /*@
97: TSTrajectoryGet - Updates the solution vector of a time stepper object by inquiring the TSTrajectory
99: Collective on TS
101: Input Parameters:
102: + tj - the trajectory object
103: . ts - the time stepper object
104: - stepnum - the step number
106: Output Parameter:
107: . time - the time associated with the step number
109: Level: developer
111: Notes: Usually one does not call this routine, it is called automatically during TSSolve()
113: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGetVecs(), TSGetSolution()
114: @*/
115: PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time)
116: {
125: if (tj->monitor) {
126: PetscViewerASCIIPrintf(tj->monitor,"TSTrajectoryGet: stepnum %D, stages %D\n",stepnum,(PetscInt)!tj->solution_only);
127: PetscViewerFlush(tj->monitor);
128: }
129: PetscLogEventBegin(TSTrajectory_Get,tj,ts,0,0);
130: (*tj->ops->get)(tj,ts,stepnum,time);
131: PetscLogEventEnd(TSTrajectory_Get,tj,ts,0,0);
132: return 0;
133: }
135: /*@
136: TSTrajectoryGetVecs - Reconstructs the vector of state and its time derivative using information from the TSTrajectory and, possibly, from the TS
138: Collective on TS
140: Input Parameters:
141: + tj - the trajectory object
142: . ts - the time stepper object (optional)
143: - stepnum - the requested step number
145: Input/Output Parameter:
147: Output Parameters:
148: + time - On input time for the step if step number is PETSC_DECIDE, on output the time associated with the step number
149: . U - state vector (can be NULL)
150: - Udot - time derivative of state vector (can be NULL)
152: Level: developer
154: Notes: If the step number is PETSC_DECIDE, the time argument is used to inquire the trajectory.
155: If the requested time does not match any in the trajectory, Lagrangian interpolations are returned.
157: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGet()
158: @*/
159: PetscErrorCode TSTrajectoryGetVecs(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time,Vec U,Vec Udot)
160: {
168: if (!U && !Udot) return 0;
170: PetscLogEventBegin(TSTrajectory_GetVecs,tj,ts,0,0);
171: if (tj->monitor) {
172: PetscInt pU,pUdot;
173: pU = U ? 1 : 0;
174: pUdot = Udot ? 1 : 0;
175: PetscViewerASCIIPrintf(tj->monitor,"Requested by GetVecs %D %D: stepnum %D, time %g\n",pU,pUdot,stepnum,(double)*time);
176: PetscViewerFlush(tj->monitor);
177: }
178: if (U && tj->lag.caching) {
179: PetscObjectId id;
180: PetscObjectState state;
182: PetscObjectStateGet((PetscObject)U,&state);
183: PetscObjectGetId((PetscObject)U,&id);
184: if (stepnum == PETSC_DECIDE) {
185: if (id == tj->lag.Ucached.id && *time == tj->lag.Ucached.time && state == tj->lag.Ucached.state) U = NULL;
186: } else {
187: if (id == tj->lag.Ucached.id && stepnum == tj->lag.Ucached.step && state == tj->lag.Ucached.state) U = NULL;
188: }
189: if (tj->monitor && !U) {
190: PetscViewerASCIIPushTab(tj->monitor);
191: PetscViewerASCIIPrintf(tj->monitor,"State vector cached\n");
192: PetscViewerASCIIPopTab(tj->monitor);
193: PetscViewerFlush(tj->monitor);
194: }
195: }
196: if (Udot && tj->lag.caching) {
197: PetscObjectId id;
198: PetscObjectState state;
200: PetscObjectStateGet((PetscObject)Udot,&state);
201: PetscObjectGetId((PetscObject)Udot,&id);
202: if (stepnum == PETSC_DECIDE) {
203: if (id == tj->lag.Udotcached.id && *time == tj->lag.Udotcached.time && state == tj->lag.Udotcached.state) Udot = NULL;
204: } else {
205: if (id == tj->lag.Udotcached.id && stepnum == tj->lag.Udotcached.step && state == tj->lag.Udotcached.state) Udot = NULL;
206: }
207: if (tj->monitor && !Udot) {
208: PetscViewerASCIIPushTab(tj->monitor);
209: PetscViewerASCIIPrintf(tj->monitor,"Derivative vector cached\n");
210: PetscViewerASCIIPopTab(tj->monitor);
211: PetscViewerFlush(tj->monitor);
212: }
213: }
214: if (!U && !Udot) {
215: PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
216: return 0;
217: }
219: if (stepnum == PETSC_DECIDE || Udot) { /* reverse search for requested time in TSHistory */
220: if (tj->monitor) {
221: PetscViewerASCIIPushTab(tj->monitor);
222: }
223: /* cached states will be updated in the function */
224: TSTrajectoryReconstruct_Private(tj,ts,*time,U,Udot);
225: if (tj->monitor) {
226: PetscViewerASCIIPopTab(tj->monitor);
227: PetscViewerFlush(tj->monitor);
228: }
229: } else if (U) { /* we were asked to load from stepnum, use TSTrajectoryGet */
230: TS fakets = ts;
231: Vec U2;
233: /* use a fake TS if ts is missing */
234: if (!ts) {
235: PetscObjectQuery((PetscObject)tj,"__fake_ts",(PetscObject*)&fakets);
236: if (!fakets) {
237: TSCreate(PetscObjectComm((PetscObject)tj),&fakets);
238: PetscObjectCompose((PetscObject)tj,"__fake_ts",(PetscObject)fakets);
239: PetscObjectDereference((PetscObject)fakets);
240: VecDuplicate(U,&U2);
241: TSSetSolution(fakets,U2);
242: PetscObjectDereference((PetscObject)U2);
243: }
244: }
245: TSTrajectoryGet(tj,fakets,stepnum,time);
246: TSGetSolution(fakets,&U2);
247: VecCopy(U2,U);
248: PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
249: PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
250: tj->lag.Ucached.time = *time;
251: tj->lag.Ucached.step = stepnum;
252: }
253: PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
254: return 0;
255: }
257: /*@C
258: TSTrajectoryViewFromOptions - View from Options
260: Collective on TSTrajectory
262: Input Parameters:
263: + A - the TSTrajectory context
264: . obj - Optional object
265: - name - command line option
267: Level: intermediate
268: .seealso: TSTrajectory, TSTrajectoryView, PetscObjectViewFromOptions(), TSTrajectoryCreate()
269: @*/
270: PetscErrorCode TSTrajectoryViewFromOptions(TSTrajectory A,PetscObject obj,const char name[])
271: {
273: PetscObjectViewFromOptions((PetscObject)A,obj,name);
274: return 0;
275: }
277: /*@C
278: TSTrajectoryView - Prints information about the trajectory object
280: Collective on TSTrajectory
282: Input Parameters:
283: + tj - the TSTrajectory context obtained from TSTrajectoryCreate()
284: - viewer - visualization context
286: Options Database Key:
287: . -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()
289: Notes:
290: The available visualization contexts include
291: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
292: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
293: output where only the first processor opens
294: the file. All other processors send their
295: data to the first processor to print.
297: The user can open an alternative visualization context with
298: PetscViewerASCIIOpen() - output to a specified file.
300: Level: developer
302: .seealso: PetscViewerASCIIOpen()
303: @*/
304: PetscErrorCode TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
305: {
306: PetscBool iascii;
309: if (!viewer) {
310: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
311: }
315: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
316: if (iascii) {
317: PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
318: PetscViewerASCIIPrintf(viewer," total number of recomputations for adjoint calculation = %D\n",tj->recomps);
319: PetscViewerASCIIPrintf(viewer," disk checkpoint reads = %D\n",tj->diskreads);
320: PetscViewerASCIIPrintf(viewer," disk checkpoint writes = %D\n",tj->diskwrites);
321: if (tj->ops->view) {
322: PetscViewerASCIIPushTab(viewer);
323: (*tj->ops->view)(tj,viewer);
324: PetscViewerASCIIPopTab(viewer);
325: }
326: }
327: return 0;
328: }
330: /*@C
331: TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory
333: Collective on TSTrajectory
335: Input Parameters:
336: + tr - the trajectory context
337: - names - the names of the components, final string must be NULL
339: Level: intermediate
341: Note: Fortran interface is not possible because of the string array argument
343: .seealso: TSTrajectory, TSGetTrajectory()
344: @*/
345: PetscErrorCode TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
346: {
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: {
372: tj->transform = transform;
373: tj->transformdestroy = destroy;
374: tj->transformctx = tctx;
375: return 0;
376: }
378: /*@
379: TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE
381: Collective
383: Input Parameter:
384: . comm - the communicator
386: Output Parameter:
387: . tj - the trajectory object
389: Level: developer
391: Notes:
392: Usually one does not call this routine, it is called automatically when one calls TSSetSaveTrajectory().
394: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
395: @*/
396: PetscErrorCode TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
397: {
398: TSTrajectory t;
401: *tj = NULL;
402: TSInitializePackage();
404: PetscHeaderCreate(t,TSTRAJECTORY_CLASSID,"TSTrajectory","Time stepping","TS",comm,TSTrajectoryDestroy,TSTrajectoryView);
405: t->setupcalled = PETSC_FALSE;
406: TSHistoryCreate(comm,&t->tsh);
408: t->lag.order = 1;
409: t->lag.L = NULL;
410: t->lag.T = NULL;
411: t->lag.W = NULL;
412: t->lag.WW = NULL;
413: t->lag.TW = NULL;
414: t->lag.TT = NULL;
415: t->lag.caching = PETSC_TRUE;
416: t->lag.Ucached.id = 0;
417: t->lag.Ucached.state = -1;
418: t->lag.Ucached.time = PETSC_MIN_REAL;
419: t->lag.Ucached.step = PETSC_MAX_INT;
420: t->lag.Udotcached.id = 0;
421: t->lag.Udotcached.state = -1;
422: t->lag.Udotcached.time = PETSC_MIN_REAL;
423: t->lag.Udotcached.step = PETSC_MAX_INT;
424: t->adjoint_solve_mode = PETSC_TRUE;
425: t->solution_only = PETSC_FALSE;
426: t->keepfiles = PETSC_FALSE;
427: t->usehistory = PETSC_TRUE;
428: *tj = t;
429: TSTrajectorySetFiletemplate(t,"TS-%06D.bin");
430: return 0;
431: }
433: /*@C
434: TSTrajectorySetType - Sets the storage method to be used as in a trajectory
436: Collective on TS
438: Input Parameters:
439: + tj - the TSTrajectory context
440: . ts - the TS context
441: - type - a known method
443: Options Database Command:
444: . -ts_trajectory_type <type> - Sets the method; use -help for a list of available methods (for instance, basic)
446: Level: developer
448: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy(), TSTrajectoryGetType()
450: @*/
451: PetscErrorCode TSTrajectorySetType(TSTrajectory tj,TS ts,TSTrajectoryType type)
452: {
453: PetscErrorCode (*r)(TSTrajectory,TS);
454: PetscBool match;
457: PetscObjectTypeCompare((PetscObject)tj,type,&match);
458: if (match) return 0;
460: PetscFunctionListFind(TSTrajectoryList,type,&r);
462: if (tj->ops->destroy) {
463: (*(tj)->ops->destroy)(tj);
465: tj->ops->destroy = NULL;
466: }
467: PetscMemzero(tj->ops,sizeof(*tj->ops));
469: PetscObjectChangeTypeName((PetscObject)tj,type);
470: (*r)(tj,ts);
471: return 0;
472: }
474: /*@C
475: TSTrajectoryGetType - Gets the trajectory type
477: Collective on TS
479: Input Parameters:
480: + tj - the TSTrajectory context
481: - ts - the TS context
483: Output Parameters:
484: . type - a known method
486: Level: developer
488: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy(), TSTrajectorySetType()
490: @*/
491: PetscErrorCode TSTrajectoryGetType(TSTrajectory tj,TS ts,TSTrajectoryType *type)
492: {
494: if (type) *type = ((PetscObject)tj)->type_name;
495: return 0;
496: }
498: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
499: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
500: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
501: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);
503: /*@C
504: TSTrajectoryRegisterAll - Registers all of the trajectory storage schecmes in the TS package.
506: Not Collective
508: Level: developer
510: .seealso: TSTrajectoryRegister()
511: @*/
512: PetscErrorCode TSTrajectoryRegisterAll(void)
513: {
514: if (TSTrajectoryRegisterAllCalled) return 0;
515: TSTrajectoryRegisterAllCalled = PETSC_TRUE;
517: TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
518: TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
519: TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
520: TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
521: return 0;
522: }
524: /*@
525: TSTrajectoryReset - Resets a trajectory context
527: Collective on TSTrajectory
529: Input Parameter:
530: . tj - the TSTrajectory context obtained from TSTrajectoryCreate()
532: Level: developer
534: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
535: @*/
536: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
537: {
538: if (!tj) return 0;
540: if (tj->ops->reset) {
541: (*tj->ops->reset)(tj);
542: }
543: PetscFree(tj->dirfiletemplate);
544: TSHistoryDestroy(&tj->tsh);
545: TSHistoryCreate(PetscObjectComm((PetscObject)tj),&tj->tsh);
546: tj->setupcalled = PETSC_FALSE;
547: return 0;
548: }
550: /*@
551: TSTrajectoryDestroy - Destroys a trajectory context
553: Collective on TSTrajectory
555: Input Parameter:
556: . tj - the TSTrajectory context obtained from TSTrajectoryCreate()
558: Level: developer
560: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
561: @*/
562: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
563: {
564: if (!*tj) return 0;
566: if (--((PetscObject)(*tj))->refct > 0) {*tj = NULL; return 0;}
568: TSTrajectoryReset(*tj);
569: TSHistoryDestroy(&(*tj)->tsh);
570: VecDestroyVecs((*tj)->lag.order+1,&(*tj)->lag.W);
571: PetscFree5((*tj)->lag.L,(*tj)->lag.T,(*tj)->lag.WW,(*tj)->lag.TT,(*tj)->lag.TW);
572: VecDestroy(&(*tj)->U);
573: VecDestroy(&(*tj)->Udot);
575: if ((*tj)->transformdestroy) (*(*tj)->transformdestroy)((*tj)->transformctx);
576: if ((*tj)->ops->destroy) (*(*tj)->ops->destroy)((*tj));
577: if (!((*tj)->keepfiles)) {
578: PetscMPIInt rank;
579: MPI_Comm comm;
581: PetscObjectGetComm((PetscObject)(*tj),&comm);
582: MPI_Comm_rank(comm,&rank);
583: if (rank == 0 && (*tj)->dirname) { /* we own the directory, so we run PetscRMTree on it */
584: PetscRMTree((*tj)->dirname);
585: }
586: }
587: PetscStrArrayDestroy(&(*tj)->names);
588: PetscFree((*tj)->dirname);
589: PetscFree((*tj)->filetemplate);
590: PetscHeaderDestroy(tj);
591: return 0;
592: }
594: /*
595: TSTrajectorySetTypeFromOptions_Private - Sets the type of ts from user options.
597: Collective on TSTrajectory
599: Input Parameter:
600: + tj - the TSTrajectory context
601: - ts - the TS context
603: Options Database Keys:
604: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
606: Level: developer
608: .seealso: TSTrajectorySetFromOptions(), TSTrajectorySetType()
609: */
610: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,TSTrajectory tj,TS ts)
611: {
612: PetscBool opt;
613: const char *defaultType;
614: char typeName[256];
616: if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
617: else defaultType = TSTRAJECTORYBASIC;
619: TSTrajectoryRegisterAll();
620: PetscOptionsFList("-ts_trajectory_type","TSTrajectory method","TSTrajectorySetType",TSTrajectoryList,defaultType,typeName,256,&opt);
621: if (opt) {
622: TSTrajectorySetType(tj,ts,typeName);
623: } else {
624: TSTrajectorySetType(tj,ts,defaultType);
625: }
626: return 0;
627: }
629: /*@
630: TSTrajectorySetUseHistory - Use TSHistory in TSTrajectory
632: Collective on TSTrajectory
634: Input Parameters:
635: + tj - the TSTrajectory context
636: - flg - PETSC_TRUE to save, PETSC_FALSE to disable
638: Options Database Keys:
639: . -ts_trajectory_use_history - have it use TSHistory
641: Level: advanced
643: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
644: @*/
645: PetscErrorCode TSTrajectorySetUseHistory(TSTrajectory tj,PetscBool flg)
646: {
649: tj->usehistory = flg;
650: return 0;
651: }
653: /*@
654: TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller
656: Collective on TSTrajectory
658: Input Parameters:
659: + tj - the TSTrajectory context
660: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
662: Options Database Keys:
663: . -ts_trajectory_monitor - print TSTrajectory information
665: Level: developer
667: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
668: @*/
669: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
670: {
673: if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
674: else tj->monitor = NULL;
675: return 0;
676: }
678: /*@
679: TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory
681: Collective on TSTrajectory
683: Input Parameters:
684: + tj - the TSTrajectory context
685: - flg - PETSC_TRUE to save, PETSC_FALSE to disable
687: Options Database Keys:
688: . -ts_trajectory_keep_files - have it keep the files
690: Notes:
691: 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.
693: Level: advanced
695: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
696: @*/
697: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
698: {
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 Parameters:
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: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
723: @*/
724: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
725: {
726: PetscBool flg;
729: PetscStrcmp(tj->dirname,dirname,&flg);
730: if (!flg && tj->dirfiletemplate) {
731: SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after TSTrajectory has been setup");
732: }
733: PetscFree(tj->dirname);
734: PetscStrallocpy(dirname,&tj->dirname);
735: return 0;
736: }
738: /*@C
739: TSTrajectorySetFiletemplate - Specify the name template for the files storing checkpoints.
741: Collective on TSTrajectory
743: Input Parameters:
744: + tj - the TSTrajectory context
745: - filetemplate - the template
747: Options Database Keys:
748: . -ts_trajectory_file_template - set the file name template
750: Notes:
751: The name template should be of the form, for example filename-%06D.bin It should not begin with a leading /
753: The final location of the files is determined by dirname/filetemplate where dirname was provided by TSTrajectorySetDirname(). The %06D is replaced by the
754: timestep counter
756: Level: developer
758: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
759: @*/
760: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
761: {
762: const char *ptr,*ptr2;
768: /* Do some cursory validation of the input. */
769: PetscStrstr(filetemplate,"%",(char**)&ptr);
771: for (ptr++; ptr && *ptr; ptr++) {
772: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
774: if (ptr2) break;
775: }
776: PetscFree(tj->filetemplate);
777: PetscStrallocpy(filetemplate,&tj->filetemplate);
778: return 0;
779: }
781: /*@
782: TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.
784: Collective on TSTrajectory
786: Input Parameters:
787: + tj - the TSTrajectory context obtained from TSTrajectoryCreate()
788: - ts - the TS context
790: Options Database Keys:
791: + -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
792: . -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
793: - -ts_trajectory_monitor - print TSTrajectory information
795: Level: developer
797: Notes:
798: This is not normally called directly by users
800: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
801: @*/
802: PetscErrorCode TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
803: {
804: PetscBool set,flg;
805: char dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];
810: PetscObjectOptionsBegin((PetscObject)tj);
811: TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
812: PetscOptionsBool("-ts_trajectory_use_history","Turn on/off usage of TSHistory",NULL,tj->usehistory,&tj->usehistory,NULL);
813: PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
814: if (set) TSTrajectorySetMonitor(tj,flg);
815: PetscOptionsInt("-ts_trajectory_reconstruction_order","Interpolation order for reconstruction",NULL,tj->lag.order,&tj->lag.order,NULL);
816: PetscOptionsBool("-ts_trajectory_reconstruction_caching","Turn on/off caching of TSTrajectoryGetVecs input",NULL,tj->lag.caching,&tj->lag.caching,NULL);
817: PetscOptionsBool("-ts_trajectory_adjointmode","Instruct the trajectory that will be used in a TSAdjointSolve()",NULL,tj->adjoint_solve_mode,&tj->adjoint_solve_mode,NULL);
818: PetscOptionsBool("-ts_trajectory_solution_only","Checkpoint solution only","TSTrajectorySetSolutionOnly",tj->solution_only,&tj->solution_only,NULL);
819: PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
820: if (set) TSTrajectorySetKeepFiles(tj,flg);
822: PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",NULL,dirname,sizeof(dirname)-14,&set);
823: if (set) {
824: TSTrajectorySetDirname(tj,dirname);
825: }
827: PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",NULL,filetemplate,sizeof(filetemplate),&set);
828: if (set) {
829: TSTrajectorySetFiletemplate(tj,filetemplate);
830: }
832: /* Handle specific TSTrajectory options */
833: if (tj->ops->setfromoptions) {
834: (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
835: }
836: PetscOptionsEnd();
837: return 0;
838: }
840: /*@
841: TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
842: of a TS trajectory.
844: Collective on TS
846: Input Parameters:
847: + ts - the TS context obtained from TSCreate()
848: - tj - the TS trajectory context
850: Level: developer
852: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
853: @*/
854: PetscErrorCode TSTrajectorySetUp(TSTrajectory tj,TS ts)
855: {
856: size_t s1,s2;
858: if (!tj) return 0;
861: if (tj->setupcalled) return 0;
863: PetscLogEventBegin(TSTrajectory_SetUp,tj,ts,0,0);
864: if (!((PetscObject)tj)->type_name) {
865: TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
866: }
867: if (tj->ops->setup) {
868: (*tj->ops->setup)(tj,ts);
869: }
871: tj->setupcalled = PETSC_TRUE;
873: /* Set the counters to zero */
874: tj->recomps = 0;
875: tj->diskreads = 0;
876: tj->diskwrites = 0;
877: PetscStrlen(tj->dirname,&s1);
878: PetscStrlen(tj->filetemplate,&s2);
879: PetscFree(tj->dirfiletemplate);
880: PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
881: PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
882: PetscLogEventEnd(TSTrajectory_SetUp,tj,ts,0,0);
883: return 0;
884: }
886: /*@
887: TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage also.
889: Collective on TSTrajectory
891: Input Parameters:
892: + tj - the TS trajectory context
893: - flg - the boolean flag
895: Level: developer
897: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
898: @*/
899: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
900: {
903: tj->solution_only = solution_only;
904: return 0;
905: }
907: /*@
908: TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.
910: Logically collective on TSTrajectory
912: Input Parameter:
913: . tj - the TS trajectory context
915: Output Parameter:
916: . flg - the boolean flag
918: Level: developer
920: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetSolutionOnly()
921: @*/
922: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj,PetscBool *solution_only)
923: {
926: *solution_only = tj->solution_only;
927: return 0;
928: }
930: /*@
931: TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.
933: Collective on TSTrajectory
935: Input Parameters:
936: + tj - the TS trajectory context
937: . ts - the TS solver context
938: - time - the requested time
940: Output Parameters:
941: + U - state vector at given time (can be interpolated)
942: - Udot - time-derivative vector at given time (can be interpolated)
944: Level: developer
946: Notes: The vectors are interpolated if time does not match any time step stored in the TSTrajectory(). Pass NULL to not request a vector.
947: This function differs from TSTrajectoryGetVecs since the vectors obtained cannot be modified, and they need to be returned by
948: calling TSTrajectoryRestoreUpdatedHistoryVecs().
950: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
951: @*/
952: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
953: {
959: if (U && !tj->U) {
960: DM dm;
962: TSGetDM(ts,&dm);
963: DMCreateGlobalVector(dm,&tj->U);
964: }
965: if (Udot && !tj->Udot) {
966: DM dm;
968: TSGetDM(ts,&dm);
969: DMCreateGlobalVector(dm,&tj->Udot);
970: }
971: TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
972: if (U) {
973: VecLockReadPush(tj->U);
974: *U = tj->U;
975: }
976: if (Udot) {
977: VecLockReadPush(tj->Udot);
978: *Udot = tj->Udot;
979: }
980: return 0;
981: }
983: /*@
984: TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with TSTrajectoryGetUpdatedHistoryVecs().
986: Collective on TSTrajectory
988: Input Parameters:
989: + tj - the TS trajectory context
990: . U - state vector at given time (can be interpolated)
991: - Udot - time-derivative vector at given time (can be interpolated)
993: Level: developer
995: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
996: @*/
997: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
998: {
1004: if (U) {
1005: VecLockReadPop(tj->U);
1006: *U = NULL;
1007: }
1008: if (Udot) {
1009: VecLockReadPop(tj->Udot);
1010: *Udot = NULL;
1011: }
1012: return 0;
1013: }