Actual source code: traj.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/tsimpl.h>
2: #include <petsc/private/tshistoryimpl.h>
3: #include <petscdm.h>
5: PetscFunctionList TSTrajectoryList = NULL;
6: PetscBool TSTrajectoryRegisterAllCalled = PETSC_FALSE;
7: PetscClassId TSTRAJECTORY_CLASSID;
8: PetscLogEvent TSTrajectory_Set, TSTrajectory_Get, TSTrajectory_GetVecs;
10: /*@C
11: TSTrajectoryRegister - Adds a way of storing trajectories to the TS package
13: Not Collective
15: Input Parameters:
16: + name - the name of a new user-defined creation routine
17: - create_func - the creation routine itself
19: Notes:
20: TSTrajectoryRegister() may be called multiple times to add several user-defined tses.
22: Level: developer
24: .keywords: TS, trajectory, timestep, register
26: .seealso: TSTrajectoryRegisterAll()
27: @*/
28: PetscErrorCode TSTrajectoryRegister(const char sname[],PetscErrorCode (*function)(TSTrajectory,TS))
29: {
33: PetscFunctionListAdd(&TSTrajectoryList,sname,function);
34: return(0);
35: }
37: /*@
38: TSTrajectorySet - Sets a vector of state in the trajectory object
40: Collective on TSTrajectory
42: Input Parameters:
43: + tj - the trajectory object
44: . ts - the time stepper object (optional)
45: . stepnum - the step number
46: . time - the current time
47: - X - the current solution
49: Level: developer
51: Notes: Usually one does not call this routine, it is called automatically during TSSolve()
53: .keywords: TS, trajectory, create
55: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectoryGet(), TSTrajectoryGetVecs()
56: @*/
57: PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
58: {
62: if (!tj) return(0);
68: if (!tj->ops->set) SETERRQ1(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"TSTrajectory type %s",((PetscObject)tj)->type_name);
69: if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
70: if (tj->monitor) {
71: PetscViewerASCIIPrintf(tj->monitor,"TSTrajectorySet: stepnum %D, time %g (stages %D)\n",stepnum,(double)time,(PetscInt)!tj->solution_only);
72: }
73: PetscLogEventBegin(TSTrajectory_Set,tj,ts,0,0);
74: (*tj->ops->set)(tj,ts,stepnum,time,X);
75: PetscLogEventEnd(TSTrajectory_Set,tj,ts,0,0);
76: TSHistoryUpdate(tj->tsh,stepnum,time);
77: if (tj->lag.caching) tj->lag.Udotcached.time = PETSC_MIN_REAL;
78: return(0);
79: }
81: /*@
82: TSTrajectoryGetNumSteps - Return the number of steps registered in the TSTrajectory via TSTrajectorySet().
84: Not collective.
86: Input Parameters:
87: . tj - the trajectory object
89: Output Parameter:
90: . steps - the number of steps
92: Level: developer
94: .keywords: TS, trajectory, create
96: .seealso: TSTrajectorySet()
97: @*/
98: PetscErrorCode TSTrajectoryGetNumSteps(TSTrajectory tj, PetscInt *steps)
99: {
105: TSHistoryGetNumSteps(tj->tsh,steps);
106: return(0);
107: }
109: /*@
110: TSTrajectoryGet - Updates the solution vector of a time stepper object by inquiring the TSTrajectory
112: Collective on TS
114: Input Parameters:
115: + tj - the trajectory object
116: . ts - the time stepper object
117: - stepnum - the step number
119: Output Parameter:
120: . time - the time associated with the step number
122: Level: developer
124: Notes: Usually one does not call this routine, it is called automatically during TSSolve()
126: .keywords: TS, trajectory, create
128: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGetVecs(), TSGetSolution()
129: @*/
130: PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time)
131: {
135: if (!tj) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS solver did not save trajectory");
140: if (!tj->ops->get) SETERRQ1(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"TSTrajectory type %s",((PetscObject)tj)->type_name);
141: if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
142: if (stepnum < 0) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"Requesting negative step number");
143: if (tj->monitor) {
144: PetscViewerASCIIPrintf(tj->monitor,"TSTrajectoryGet: stepnum %D, stages %D\n",stepnum,(PetscInt)!tj->solution_only);
145: PetscViewerFlush(tj->monitor);
146: }
147: PetscLogEventBegin(TSTrajectory_Get,tj,ts,0,0);
148: (*tj->ops->get)(tj,ts,stepnum,time);
149: PetscLogEventEnd(TSTrajectory_Get,tj,ts,0,0);
150: return(0);
151: }
153: /*@
154: TSTrajectoryGetVecs - Reconstructs the vector of state and its time derivative using information from the TSTrajectory and, possibly, from the TS
156: Collective on TS
158: Input Parameters:
159: + tj - the trajectory object
160: . ts - the time stepper object (optional)
161: - stepnum - the requested step number
163: Input/Output Parameters:
164: . time - the time associated with the step number
166: Output Parameters:
167: + U - state vector (can be NULL)
168: - Udot - time derivative of state vector (can be NULL)
170: Level: developer
172: Notes: If the step number is PETSC_DECIDE, the time argument is used to inquire the trajectory.
173: If the requested time does not match any in the trajectory, Lagrangian interpolations are returned.
175: .keywords: TS, trajectory, create
177: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySet(), TSTrajectoryGet()
178: @*/
179: PetscErrorCode TSTrajectoryGetVecs(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time,Vec U,Vec Udot)
180: {
184: if (!tj) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS solver did not save trajectory");
191: if (!U && !Udot) return(0);
192: if (!tj->setupcalled) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ORDER,"TSTrajectorySetUp should be called first");
193: PetscLogEventBegin(TSTrajectory_GetVecs,tj,ts,0,0);
194: if (tj->monitor) {
195: PetscInt pU,pUdot;
196: pU = U ? 1 : 0;
197: pUdot = Udot ? 1 : 0;
198: PetscViewerASCIIPrintf(tj->monitor,"Requested by GetVecs %D %D: stepnum %D, time %g\n",pU,pUdot,stepnum,(double)*time);
199: PetscViewerFlush(tj->monitor);
200: }
201: if (U && tj->lag.caching) {
202: PetscObjectId id;
203: PetscObjectState state;
205: PetscObjectStateGet((PetscObject)U,&state);
206: PetscObjectGetId((PetscObject)U,&id);
207: if (stepnum == PETSC_DECIDE) {
208: if (id == tj->lag.Ucached.id && *time == tj->lag.Ucached.time && state == tj->lag.Ucached.state) U = NULL;
209: } else {
210: if (id == tj->lag.Ucached.id && stepnum == tj->lag.Ucached.step && state == tj->lag.Ucached.state) U = NULL;
211: }
212: if (tj->monitor && !U) {
213: PetscViewerASCIIPushTab(tj->monitor);
214: PetscViewerASCIIPrintf(tj->monitor,"State vector cached\n");
215: PetscViewerASCIIPopTab(tj->monitor);
216: PetscViewerFlush(tj->monitor);
217: }
218: }
219: if (Udot && tj->lag.caching) {
220: PetscObjectId id;
221: PetscObjectState state;
223: PetscObjectStateGet((PetscObject)Udot,&state);
224: PetscObjectGetId((PetscObject)Udot,&id);
225: if (stepnum == PETSC_DECIDE) {
226: if (id == tj->lag.Udotcached.id && *time == tj->lag.Udotcached.time && state == tj->lag.Udotcached.state) Udot = NULL;
227: } else {
228: if (id == tj->lag.Udotcached.id && stepnum == tj->lag.Udotcached.step && state == tj->lag.Udotcached.state) Udot = NULL;
229: }
230: if (tj->monitor && !Udot) {
231: PetscViewerASCIIPushTab(tj->monitor);
232: PetscViewerASCIIPrintf(tj->monitor,"Derivative vector cached\n");
233: PetscViewerASCIIPopTab(tj->monitor);
234: PetscViewerFlush(tj->monitor);
235: }
236: }
237: if (!U && !Udot) {
238: PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
239: return(0);
240: }
242: if (stepnum == PETSC_DECIDE || Udot) { /* reverse search for requested time in TSHistory */
243: if (tj->monitor) {
244: PetscViewerASCIIPushTab(tj->monitor);
245: }
246: /* cached states will be updated in the function */
247: TSTrajectoryReconstruct_Private(tj,ts,*time,U,Udot);
248: if (tj->monitor) {
249: PetscViewerASCIIPopTab(tj->monitor);
250: PetscViewerFlush(tj->monitor);
251: }
252: } else if (U) { /* we were asked to load from stepnum, use TSTrajectoryGet */
253: TS fakets = ts;
254: Vec U2;
256: /* use a fake TS if ts is missing */
257: if (!ts) {
258: PetscObjectQuery((PetscObject)tj,"__fake_ts",(PetscObject*)&fakets);
259: if (!fakets) {
260: TSCreate(PetscObjectComm((PetscObject)tj),&fakets);
261: PetscObjectCompose((PetscObject)tj,"__fake_ts",(PetscObject)fakets);
262: PetscObjectDereference((PetscObject)fakets);
263: VecDuplicate(U,&U2);
264: TSSetSolution(fakets,U2);
265: PetscObjectDereference((PetscObject)U2);
266: }
267: }
268: TSTrajectoryGet(tj,fakets,stepnum,time);
269: TSGetSolution(fakets,&U2);
270: VecCopy(U2,U);
271: PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
272: PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
273: tj->lag.Ucached.time = *time;
274: tj->lag.Ucached.step = stepnum;
275: }
276: PetscLogEventEnd(TSTrajectory_GetVecs,tj,ts,0,0);
277: return(0);
278: }
280: /*@C
281: TSTrajectoryView - Prints information about the trajectory object
283: Collective on TSTrajectory
285: Input Parameters:
286: + tj - the TSTrajectory context obtained from TSTrajectoryCreate()
287: - viewer - visualization context
289: Options Database Key:
290: . -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()
292: Notes:
293: The available visualization contexts include
294: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
295: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
296: output where only the first processor opens
297: the file. All other processors send their
298: data to the first processor to print.
300: The user can open an alternative visualization context with
301: PetscViewerASCIIOpen() - output to a specified file.
303: Level: developer
305: .keywords: TS, trajectory, timestep, view
307: .seealso: PetscViewerASCIIOpen()
308: @*/
309: PetscErrorCode TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
310: {
312: PetscBool iascii;
316: if (!viewer) {
317: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
318: }
322: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
323: if (iascii) {
324: PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
325: PetscViewerASCIIPrintf(viewer," total number of recomputations for adjoint calculation = %D\n",tj->recomps);
326: PetscViewerASCIIPrintf(viewer," disk checkpoint reads = %D\n",tj->diskreads);
327: PetscViewerASCIIPrintf(viewer," disk checkpoint writes = %D\n",tj->diskwrites);
328: if (tj->ops->view) {
329: PetscViewerASCIIPushTab(viewer);
330: (*tj->ops->view)(tj,viewer);
331: PetscViewerASCIIPopTab(viewer);
332: }
333: }
334: return(0);
335: }
337: /*@C
338: TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory
340: Collective on TSTrajectory
342: Input Parameters:
343: + tr - the trajectory context
344: - names - the names of the components, final string must be NULL
346: Level: intermediate
348: Note: Fortran interface is not possible because of the string array argument
350: .keywords: TS, TSTrajectory, vector, monitor, view
352: .seealso: TSTrajectory, TSGetTrajectory()
353: @*/
354: PetscErrorCode TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
355: {
356: PetscErrorCode ierr;
361: PetscStrArrayDestroy(&ctx->names);
362: PetscStrArrayallocpy(names,&ctx->names);
363: return(0);
364: }
366: /*@C
367: TSTrajectorySetTransform - Solution vector will be transformed by provided function before being saved to disk
369: Collective on TSLGCtx
371: Input Parameters:
372: + tj - the TSTrajectory context
373: . transform - the transform function
374: . destroy - function to destroy the optional context
375: - ctx - optional context used by transform function
377: Level: intermediate
379: .keywords: TSTrajectory, vector, monitor, view
381: .seealso: TSTrajectorySetVariableNames(), TSTrajectory, TSMonitorLGSetTransform()
382: @*/
383: PetscErrorCode TSTrajectorySetTransform(TSTrajectory tj,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
384: {
387: tj->transform = transform;
388: tj->transformdestroy = destroy;
389: tj->transformctx = tctx;
390: return(0);
391: }
393: /*@
394: TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE
396: Collective on MPI_Comm
398: Input Parameter:
399: . comm - the communicator
401: Output Parameter:
402: . tj - the trajectory object
404: Level: developer
406: Notes:
407: Usually one does not call this routine, it is called automatically when one calls TSSetSaveTrajectory().
409: .keywords: TS, trajectory, create
411: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
412: @*/
413: PetscErrorCode TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
414: {
415: TSTrajectory t;
420: *tj = NULL;
421: TSInitializePackage();
423: PetscHeaderCreate(t,TSTRAJECTORY_CLASSID,"TSTrajectory","Time stepping","TS",comm,TSTrajectoryDestroy,TSTrajectoryView);
424: t->setupcalled = PETSC_FALSE;
425: TSHistoryCreate(comm,&t->tsh);
427: t->lag.order = 1;
428: t->lag.L = NULL;
429: t->lag.T = NULL;
430: t->lag.W = NULL;
431: t->lag.WW = NULL;
432: t->lag.TW = NULL;
433: t->lag.TT = NULL;
434: t->lag.caching = PETSC_TRUE;
435: t->lag.Ucached.id = 0;
436: t->lag.Ucached.state = -1;
437: t->lag.Ucached.time = PETSC_MIN_REAL;
438: t->lag.Ucached.step = PETSC_MAX_INT;
439: t->lag.Udotcached.id = 0;
440: t->lag.Udotcached.state = -1;
441: t->lag.Udotcached.time = PETSC_MIN_REAL;
442: t->lag.Udotcached.step = PETSC_MAX_INT;
443: t->adjoint_solve_mode = PETSC_TRUE;
444: t->solution_only = PETSC_FALSE;
445: t->keepfiles = PETSC_FALSE;
446: *tj = t;
447: TSTrajectorySetDirname(t,"SA-data");
448: TSTrajectorySetFiletemplate(t,"SA-%06D.bin");
449: return(0);
450: }
452: /*@C
453: TSTrajectorySetType - Sets the storage method to be used as in a trajectory
455: Collective on TS
457: Input Parameters:
458: + tj - the TSTrajectory context
459: . ts - the TS context
460: - type - a known method
462: Options Database Command:
463: . -ts_trajectory_type <type> - Sets the method; use -help for a list of available methods (for instance, basic)
465: Level: developer
467: .keywords: TS, trajectory, timestep, set, type
469: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy()
471: @*/
472: PetscErrorCode TSTrajectorySetType(TSTrajectory tj,TS ts,TSTrajectoryType type)
473: {
474: PetscErrorCode (*r)(TSTrajectory,TS);
475: PetscBool match;
480: PetscObjectTypeCompare((PetscObject)tj,type,&match);
481: if (match) return(0);
483: PetscFunctionListFind(TSTrajectoryList,type,&r);
484: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSTrajectory type: %s",type);
485: if (tj->ops->destroy) {
486: (*(tj)->ops->destroy)(tj);
488: tj->ops->destroy = NULL;
489: }
490: PetscMemzero(tj->ops,sizeof(*tj->ops));
492: PetscObjectChangeTypeName((PetscObject)tj,type);
493: (*r)(tj,ts);
494: return(0);
495: }
497: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
498: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
499: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
500: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);
502: /*@C
503: TSTrajectoryRegisterAll - Registers all of the trajectory storage schecmes in the TS package.
505: Not Collective
507: Level: developer
509: .keywords: TS, trajectory, register, all
511: .seealso: TSTrajectoryRegister()
512: @*/
513: PetscErrorCode TSTrajectoryRegisterAll(void)
514: {
518: if (TSTrajectoryRegisterAllCalled) return(0);
519: TSTrajectoryRegisterAllCalled = PETSC_TRUE;
521: TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
522: TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
523: TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
524: TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
525: return(0);
526: }
528: /*@
529: TSTrajectoryReset - Resets a trajectory context
531: Collective on TSTrajectory
533: Input Parameter:
534: . tj - the TSTrajectory context obtained from TSTrajectoryCreate()
536: Level: developer
538: .keywords: TS, trajectory, timestep, reset
540: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
541: @*/
542: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
543: {
547: if (!tj) return(0);
549: if (tj->ops->reset) {
550: (*tj->ops->reset)(tj);
551: }
552: PetscFree(tj->dirfiletemplate);
553: TSHistoryDestroy(&tj->tsh);
554: TSHistoryCreate(PetscObjectComm((PetscObject)tj),&tj->tsh);
555: tj->setupcalled = PETSC_FALSE;
556: return(0);
557: }
559: /*@
560: TSTrajectoryDestroy - Destroys a trajectory context
562: Collective on TSTrajectory
564: Input Parameter:
565: . tj - the TSTrajectory context obtained from TSTrajectoryCreate()
567: Level: developer
569: .keywords: TS, trajectory, timestep, destroy
571: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
572: @*/
573: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
574: {
578: if (!*tj) return(0);
580: if (--((PetscObject)(*tj))->refct > 0) {*tj = 0; return(0);}
582: TSTrajectoryReset(*tj);
583: TSHistoryDestroy(&(*tj)->tsh);
584: VecDestroyVecs((*tj)->lag.order+1,&(*tj)->lag.W);
585: PetscFree5((*tj)->lag.L,(*tj)->lag.T,(*tj)->lag.WW,(*tj)->lag.TT,(*tj)->lag.TW);
586: VecDestroy(&(*tj)->U);
587: VecDestroy(&(*tj)->Udot);
589: if ((*tj)->transformdestroy) {(*(*tj)->transformdestroy)((*tj)->transformctx);}
590: if ((*tj)->ops->destroy) {(*(*tj)->ops->destroy)((*tj));}
591: if (!((*tj)->keepfiles)) {
592: PetscMPIInt rank;
593: MPI_Comm comm;
595: PetscObjectGetComm((PetscObject)(*tj),&comm);
596: MPI_Comm_rank(comm,&rank);
597: if (!rank) { /* we own the directory, so we run PetscRMTree on it */
598: PetscRMTree((*tj)->dirname);
599: }
600: }
601: PetscStrArrayDestroy(&(*tj)->names);
602: PetscFree((*tj)->dirname);
603: PetscFree((*tj)->filetemplate);
604: PetscHeaderDestroy(tj);
605: return(0);
606: }
608: /*
609: TSTrajectorySetTypeFromOptions_Private - Sets the type of ts from user options.
611: Collective on TSTrajectory
613: Input Parameter:
614: + tj - the TSTrajectory context
615: - ts - the TS context
617: Options Database Keys:
618: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
620: Level: developer
622: .keywords: TS, trajectory, set, options, type
624: .seealso: TSTrajectorySetFromOptions(), TSTrajectorySetType()
625: */
626: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,TSTrajectory tj,TS ts)
627: {
628: PetscBool opt;
629: const char *defaultType;
630: char typeName[256];
634: if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
635: else defaultType = TSTRAJECTORYBASIC;
637: TSTrajectoryRegisterAll();
638: PetscOptionsFList("-ts_trajectory_type","TSTrajectory method","TSTrajectorySetType",TSTrajectoryList,defaultType,typeName,256,&opt);
639: if (opt) {
640: TSTrajectorySetType(tj,ts,typeName);
641: } else {
642: TSTrajectorySetType(tj,ts,defaultType);
643: }
644: return(0);
645: }
647: /*@
648: TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller
650: Collective on TSTrajectory
652: Input Arguments:
653: + tj - the TSTrajectory context
654: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
656: Options Database Keys:
657: . -ts_trajectory_monitor - print TSTrajectory information
659: Level: developer
661: .keywords: TS, trajectory, set, monitor
663: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
664: @*/
665: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
666: {
670: if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
671: else tj->monitor = NULL;
672: return(0);
673: }
675: /*@
676: TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory
678: Collective on TSTrajectory
680: Input Arguments:
681: + tj - the TSTrajectory context
682: - flg - PETSC_TRUE to save, PETSC_FALSE to disable
684: Options Database Keys:
685: . -ts_trajectory_keep_files - have it keep the files
687: Notes:
688: By default the TSTrajectory used for adjoint computations, TSTRAJECTORYBASIC, removes the files it generates at the end of the run. This causes the files to be kept.
690: Level: advanced
692: .keywords: TS, trajectory, set, monitor
694: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
695: @*/
696: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
697: {
701: tj->keepfiles = flg;
702: return(0);
703: }
705: /*@C
706: TSTrajectorySetDirname - Specify the name of the directory where disk checkpoints are stored.
708: Collective on TSTrajectory
710: Input Arguments:
711: + tj - the TSTrajectory context
712: - dirname - the directory name
714: Options Database Keys:
715: . -ts_trajectory_dirname - set the directory name
717: Notes:
718: The final location of the files is determined by dirname/filetemplate where filetemplate was provided by TSTrajectorySetFiletemplate()
720: Level: developer
722: .keywords: TS, trajectory, set
724: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
725: @*/
726: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
727: {
729: PetscBool flg;
733: PetscStrcmp(tj->dirname,dirname,&flg);
734: if (!flg && tj->dirfiletemplate) {
735: SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after TSTrajectory has been setup");
736: }
737: PetscFree(tj->dirname);
738: PetscStrallocpy(dirname,&tj->dirname);
739: return(0);
740: }
742: /*@C
743: TSTrajectorySetFiletemplate - Specify the name template for the files storing checkpoints.
745: Collective on TSTrajectory
747: Input Arguments:
748: + tj - the TSTrajectory context
749: - filetemplate - the template
751: Options Database Keys:
752: . -ts_trajectory_file_template - set the file name template
754: Notes:
755: The name template should be of the form, for example filename-%06D.bin It should not begin with a leading /
757: The final location of the files is determined by dirname/filetemplate where dirname was provided by TSTrajectorySetDirname(). The %06D is replaced by the
758: timestep counter
760: Level: developer
762: .keywords: TS, trajectory, set
764: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
765: @*/
766: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
767: {
769: const char *ptr,*ptr2;
773: if (tj->dirfiletemplate) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set filetemplate after TSTrajectory has been setup");
775: if (!filetemplate[0]) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
776: /* Do some cursory validation of the input. */
777: PetscStrstr(filetemplate,"%",(char**)&ptr);
778: if (!ptr) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
779: for (ptr++; ptr && *ptr; ptr++) {
780: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
781: if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"Invalid file template argument to -ts_trajectory_file_template, should look like filename-%%06D.bin");
782: if (ptr2) break;
783: }
784: PetscFree(tj->filetemplate);
785: PetscStrallocpy(filetemplate,&tj->filetemplate);
786: return(0);
787: }
789: /*@
790: TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.
792: Collective on TSTrajectory
794: Input Parameter:
795: + tj - the TSTrajectory context obtained from TSTrajectoryCreate()
796: - ts - the TS context
798: Options Database Keys:
799: + -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
800: . -ts_trajectory_keep_files <true,false> - keep the files generated by the code after the program ends. This is true by default for TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
801: - -ts_trajectory_monitor - print TSTrajectory information
803: Level: developer
805: Notes:
806: This is not normally called directly by users
808: .keywords: TS, trajectory, timestep, set, options, database
810: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
811: @*/
812: PetscErrorCode TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
813: {
814: PetscBool set,flg;
815: char dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];
821: PetscObjectOptionsBegin((PetscObject)tj);
822: TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
823: PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
824: if (set) {TSTrajectorySetMonitor(tj,flg);}
825: PetscOptionsInt("-ts_trajectory_reconstruction_order","Interpolation order for reconstruction",NULL,tj->lag.order,&tj->lag.order,NULL);
826: PetscOptionsBool("-ts_trajectory_reconstruction_caching","Turn on/off caching of TSTrajectoryGetVecs input",NULL,tj->lag.caching,&tj->lag.caching,NULL);
827: PetscOptionsBool("-ts_trajectory_adjointmode","Instruct the trajectory that will be used in a TSAdjointSolve()",NULL,tj->adjoint_solve_mode,&tj->adjoint_solve_mode,NULL);
828: PetscOptionsBool("-ts_trajectory_solution_only","Checkpoint solution only","TSTrajectorySetSolutionOnly",tj->solution_only,&tj->solution_only,NULL);
829: PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
830: if (set) {TSTrajectorySetKeepFiles(tj,flg);}
832: PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",0,dirname,PETSC_MAX_PATH_LEN-14,&set);
833: if (set) {
834: TSTrajectorySetDirname(tj,dirname);
835: }
837: PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",0,filetemplate,PETSC_MAX_PATH_LEN,&set);
838: if (set) {
839: TSTrajectorySetFiletemplate(tj,filetemplate);
840: }
842: /* Handle specific TSTrajectory options */
843: if (tj->ops->setfromoptions) {
844: (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
845: }
846: PetscOptionsEnd();
847: return(0);
848: }
850: /*@
851: TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
852: of a TS trajectory.
854: Collective on TS
856: Input Parameter:
857: + ts - the TS context obtained from TSCreate()
858: - tj - the TS trajectory context
860: Level: developer
862: .keywords: TS, trajectory, setup
864: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
865: @*/
866: PetscErrorCode TSTrajectorySetUp(TSTrajectory tj,TS ts)
867: {
869: size_t s1,s2;
872: if (!tj) return(0);
875: if (tj->setupcalled) return(0);
877: if (!((PetscObject)tj)->type_name) {
878: TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
879: }
880: if (tj->ops->setup) {
881: (*tj->ops->setup)(tj,ts);
882: }
884: tj->setupcalled = PETSC_TRUE;
886: /* Set the counters to zero */
887: tj->recomps = 0;
888: tj->diskreads = 0;
889: tj->diskwrites = 0;
890: PetscStrlen(tj->dirname,&s1);
891: PetscStrlen(tj->filetemplate,&s2);
892: PetscFree(tj->dirfiletemplate);
893: PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
894: PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
895: return(0);
896: }
898: /*@
899: TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage also.
901: Collective on TSTrajectory
903: Input Parameter:
904: + tj - the TS trajectory context
905: - flg - the boolean flag
907: Level: developer
909: .keywords: trajectory
911: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryGetSolutionOnly()
912: @*/
913: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj,PetscBool solution_only)
914: {
918: tj->solution_only = solution_only;
919: return(0);
920: }
922: /*@
923: TSTrajectoryGetSolutionOnly - Gets the value set with TSTrajectorySetSolutionOnly.
925: Logically collective on TSTrajectory
927: Input Parameter:
928: . tj - the TS trajectory context
930: Output Parameter:
931: - flg - the boolean flag
933: Level: developer
935: .keywords: trajectory
937: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetSolutionOnly()
938: @*/
939: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj,PetscBool *solution_only)
940: {
944: *solution_only = tj->solution_only;
945: return(0);
946: }
948: /*@
949: TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.
951: Collective on TSTrajectory
953: Input Parameter:
954: + tj - the TS trajectory context
955: . ts - the TS solver context
956: - time - the requested time
958: Output Parameter:
959: + U - state vector at given time (can be interpolated)
960: - Udot - time-derivative vector at given time (can be interpolated)
962: Level: developer
964: Notes: The vectors are interpolated if time does not match any time step stored in the TSTrajectory(). Pass NULL to not request a vector.
965: This function differs from TSTrajectoryGetVecs since the vectors obtained cannot be modified, and they need to be returned by
966: calling TSTrajectoryRestoreUpdatedHistoryVecs().
968: .keywords: trajectory
970: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectoryRestoreUpdatedHistoryVecs(), TSTrajectoryGetVecs()
971: @*/
972: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
973: {
982: if (U && !tj->U) {
983: DM dm;
985: TSGetDM(ts,&dm);
986: DMCreateGlobalVector(dm,&tj->U);
987: }
988: if (Udot && !tj->Udot) {
989: DM dm;
991: TSGetDM(ts,&dm);
992: DMCreateGlobalVector(dm,&tj->Udot);
993: }
994: TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&time,U ? tj->U : NULL,Udot ? tj->Udot : NULL);
995: if (U) {
996: VecLockReadPush(tj->U);
997: *U = tj->U;
998: }
999: if (Udot) {
1000: VecLockReadPush(tj->Udot);
1001: *Udot = tj->Udot;
1002: }
1003: return(0);
1004: }
1006: /*@
1007: TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with TSTrajectoryGetUpdatedHistoryVecs().
1009: Collective on TSTrajectory
1011: Input Parameter:
1012: + tj - the TS trajectory context
1013: . U - state vector at given time (can be interpolated)
1014: - Udot - time-derivative vector at given time (can be interpolated)
1016: Level: developer
1018: .keywords: trajectory
1020: .seealso: TSTrajectoryGetUpdatedHistoryVecs()
1021: @*/
1022: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1023: {
1030: if (U && *U != tj->U) SETERRQ(PetscObjectComm((PetscObject)*U),PETSC_ERR_USER,"U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1031: if (Udot && *Udot != tj->Udot) SETERRQ(PetscObjectComm((PetscObject)*Udot),PETSC_ERR_USER,"Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1032: if (U) {
1033: VecLockReadPop(tj->U);
1034: *U = NULL;
1035: }
1036: if (Udot) {
1037: VecLockReadPop(tj->Udot);
1038: *Udot = NULL;
1039: }
1040: return(0);
1041: }