Actual source code: traj.c
petsc-3.12.5 2020-03-29
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: }