Actual source code: tsmon.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdm.h>
3: #include <petscds.h>
4: #include <petscdmswarm.h>
5: #include <petscdraw.h>
7: /*@C
8: TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()
10: Collective on TS
12: Input Parameters:
13: + ts - time stepping context obtained from TSCreate()
14: . step - step number that has just completed
15: . ptime - model time of the state
16: - u - state at the current model time
18: Notes:
19: TSMonitor() is typically used automatically within the time stepping implementations.
20: Users would almost never call this routine directly.
22: A step of -1 indicates that the monitor is being called on a solution obtained by interpolating from computed solutions
24: Level: developer
26: @*/
27: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
28: {
29: DM dm;
30: PetscInt i,n = ts->numbermonitors;
35: TSGetDM(ts,&dm);
36: DMSetOutputSequenceNumber(dm,step,ptime);
38: VecLockReadPush(u);
39: for (i=0; i<n; i++) {
40: (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
41: }
42: VecLockReadPop(u);
43: return 0;
44: }
46: /*@C
47: TSMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
49: Collective on TS
51: Input Parameters:
52: + ts - TS object you wish to monitor
53: . name - the monitor type one is seeking
54: . help - message indicating what monitoring is done
55: . manual - manual page for the monitor
56: . monitor - the monitor function
57: - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the TS or PetscViewer objects
59: Level: developer
61: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
62: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
63: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
64: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
65: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
66: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
67: PetscOptionsFList(), PetscOptionsEList()
68: @*/
69: PetscErrorCode TSMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
70: {
71: PetscViewer viewer;
72: PetscViewerFormat format;
73: PetscBool flg;
75: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
76: if (flg) {
77: PetscViewerAndFormat *vf;
78: PetscViewerAndFormatCreate(viewer,format,&vf);
79: PetscObjectDereference((PetscObject)viewer);
80: if (monitorsetup) {
81: (*monitorsetup)(ts,vf);
82: }
83: TSMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
84: }
85: return 0;
86: }
88: /*@C
89: TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
90: timestep to display the iteration's progress.
92: Logically Collective on TS
94: Input Parameters:
95: + ts - the TS context obtained from TSCreate()
96: . monitor - monitoring routine
97: . mctx - [optional] user-defined context for private data for the
98: monitor routine (use NULL if no context is desired)
99: - monitordestroy - [optional] routine that frees monitor context
100: (may be NULL)
102: Calling sequence of monitor:
103: $ PetscErrorCode monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
105: + ts - the TS context
106: . steps - iteration number (after the final time step the monitor routine may be called with a step of -1, this indicates the solution has been interpolated to this time)
107: . time - current time
108: . u - current iterate
109: - mctx - [optional] monitoring context
111: Notes:
112: This routine adds an additional monitor to the list of monitors that
113: already has been loaded.
115: Fortran Notes:
116: Only a single monitor function can be set for each TS object
118: Level: intermediate
120: .seealso: TSMonitorDefault(), TSMonitorCancel()
121: @*/
122: PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
123: {
124: PetscInt i;
125: PetscBool identical;
128: for (i=0; i<ts->numbermonitors;i++) {
129: PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,mdestroy,(PetscErrorCode (*)(void))ts->monitor[i],ts->monitorcontext[i],ts->monitordestroy[i],&identical);
130: if (identical) return 0;
131: }
133: ts->monitor[ts->numbermonitors] = monitor;
134: ts->monitordestroy[ts->numbermonitors] = mdestroy;
135: ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
136: return 0;
137: }
139: /*@C
140: TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
142: Logically Collective on TS
144: Input Parameters:
145: . ts - the TS context obtained from TSCreate()
147: Notes:
148: There is no way to remove a single, specific monitor.
150: Level: intermediate
152: .seealso: TSMonitorDefault(), TSMonitorSet()
153: @*/
154: PetscErrorCode TSMonitorCancel(TS ts)
155: {
156: PetscInt i;
159: for (i=0; i<ts->numbermonitors; i++) {
160: if (ts->monitordestroy[i]) {
161: (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
162: }
163: }
164: ts->numbermonitors = 0;
165: return 0;
166: }
168: /*@C
169: TSMonitorDefault - The Default monitor, prints the timestep and time for each step
171: Level: intermediate
173: .seealso: TSMonitorSet()
174: @*/
175: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
176: {
177: PetscViewer viewer = vf->viewer;
178: PetscBool iascii,ibinary;
181: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
182: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
183: PetscViewerPushFormat(viewer,vf->format);
184: if (iascii) {
185: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
186: if (step == -1) { /* this indicates it is an interpolated solution */
187: PetscViewerASCIIPrintf(viewer,"Interpolated solution at time %g between steps %D and %D\n",(double)ptime,ts->steps-1,ts->steps);
188: } else {
189: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
190: }
191: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
192: } else if (ibinary) {
193: PetscMPIInt rank;
194: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);
195: if (!rank) {
196: PetscBool skipHeader;
197: PetscInt classid = REAL_FILE_CLASSID;
199: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
200: if (!skipHeader) {
201: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
202: }
203: PetscRealView(1,&ptime,viewer);
204: } else {
205: PetscRealView(0,&ptime,viewer);
206: }
207: }
208: PetscViewerPopFormat(viewer);
209: return 0;
210: }
212: /*@C
213: TSMonitorExtreme - Prints the extreme values of the solution at each timestep
215: Level: intermediate
217: .seealso: TSMonitorSet()
218: @*/
219: PetscErrorCode TSMonitorExtreme(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
220: {
221: PetscViewer viewer = vf->viewer;
222: PetscBool iascii;
223: PetscReal max,min;
226: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
227: PetscViewerPushFormat(viewer,vf->format);
228: if (iascii) {
229: VecMax(v,NULL,&max);
230: VecMin(v,NULL,&min);
231: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
232: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s max %g min %g\n",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)" : "",(double)max,(double)min);
233: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
234: }
235: PetscViewerPopFormat(viewer);
236: return 0;
237: }
239: /*@C
240: TSMonitorLGCtxCreate - Creates a TSMonitorLGCtx context for use with
241: TS to monitor the solution process graphically in various ways
243: Collective on TS
245: Input Parameters:
246: + host - the X display to open, or null for the local machine
247: . label - the title to put in the title bar
248: . x, y - the screen coordinates of the upper left coordinate of the window
249: . m, n - the screen width and height in pixels
250: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
252: Output Parameter:
253: . ctx - the context
255: Options Database Key:
256: + -ts_monitor_lg_timestep - automatically sets line graph monitor
257: + -ts_monitor_lg_timestep_log - automatically sets line graph monitor
258: . -ts_monitor_lg_solution - monitor the solution (or certain values of the solution by calling TSMonitorLGSetDisplayVariables() or TSMonitorLGCtxSetDisplayVariables())
259: . -ts_monitor_lg_error - monitor the error
260: . -ts_monitor_lg_ksp_iterations - monitor the number of KSP iterations needed for each timestep
261: . -ts_monitor_lg_snes_iterations - monitor the number of SNES iterations needed for each timestep
262: - -lg_use_markers <true,false> - mark the data points (at each time step) on the plot; default is true
264: Notes:
265: Use TSMonitorLGCtxDestroy() to destroy.
267: One can provide a function that transforms the solution before plotting it with TSMonitorLGCtxSetTransform() or TSMonitorLGSetTransform()
269: Many of the functions that control the monitoring have two forms: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a TS object as the
270: first argument (if that TS object does not have a TSMonitorLGCtx associated with it the function call is ignored) and the second takes a TSMonitorLGCtx object
271: as the first argument.
273: One can control the names displayed for each solution or error variable with TSMonitorLGCtxSetVariableNames() or TSMonitorLGSetVariableNames()
275: Level: intermediate
277: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError(), TSMonitorDefault(), VecView(),
278: TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
279: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
280: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
281: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
283: @*/
284: PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
285: {
286: PetscDraw draw;
288: PetscNew(ctx);
289: PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
290: PetscDrawSetFromOptions(draw);
291: PetscDrawLGCreate(draw,1,&(*ctx)->lg);
292: PetscDrawLGSetFromOptions((*ctx)->lg);
293: PetscDrawDestroy(&draw);
294: (*ctx)->howoften = howoften;
295: return 0;
296: }
298: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
299: {
300: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
301: PetscReal x = ptime,y;
303: if (step < 0) return 0; /* -1 indicates an interpolated solution */
304: if (!step) {
305: PetscDrawAxis axis;
306: const char *ylabel = ctx->semilogy ? "Log Time Step" : "Time Step";
307: PetscDrawLGGetAxis(ctx->lg,&axis);
308: PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time",ylabel);
309: PetscDrawLGReset(ctx->lg);
310: }
311: TSGetTimeStep(ts,&y);
312: if (ctx->semilogy) y = PetscLog10Real(y);
313: PetscDrawLGAddPoint(ctx->lg,&x,&y);
314: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
315: PetscDrawLGDraw(ctx->lg);
316: PetscDrawLGSave(ctx->lg);
317: }
318: return 0;
319: }
321: /*@C
322: TSMonitorLGCtxDestroy - Destroys a line graph context that was created
323: with TSMonitorLGCtxCreate().
325: Collective on TSMonitorLGCtx
327: Input Parameter:
328: . ctx - the monitor context
330: Level: intermediate
332: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep();
333: @*/
334: PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
335: {
336: if ((*ctx)->transformdestroy) {
337: ((*ctx)->transformdestroy)((*ctx)->transformctx);
338: }
339: PetscDrawLGDestroy(&(*ctx)->lg);
340: PetscStrArrayDestroy(&(*ctx)->names);
341: PetscStrArrayDestroy(&(*ctx)->displaynames);
342: PetscFree((*ctx)->displayvariables);
343: PetscFree((*ctx)->displayvalues);
344: PetscFree(*ctx);
345: return 0;
346: }
348: /* Creates a TS Monitor SPCtx for use with DMSwarm particle visualizations */
349: PetscErrorCode TSMonitorSPCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,PetscInt retain,PetscBool phase,TSMonitorSPCtx *ctx)
350: {
351: PetscDraw draw;
353: PetscNew(ctx);
354: PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
355: PetscDrawSetFromOptions(draw);
356: PetscDrawSPCreate(draw,1,&(*ctx)->sp);
357: PetscDrawDestroy(&draw);
358: (*ctx)->howoften = howoften;
359: (*ctx)->retain = retain;
360: (*ctx)->phase = phase;
361: return 0;
362: }
364: /*
365: Destroys a TSMonitorSPCtx that was created with TSMonitorSPCtxCreate
366: */
367: PetscErrorCode TSMonitorSPCtxDestroy(TSMonitorSPCtx *ctx)
368: {
370: PetscDrawSPDestroy(&(*ctx)->sp);
371: PetscFree(*ctx);
373: return 0;
375: }
377: /*@C
378: TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
379: VecView() for the solution at each timestep
381: Collective on TS
383: Input Parameters:
384: + ts - the TS context
385: . step - current time-step
386: . ptime - current time
387: - dummy - either a viewer or NULL
389: Options Database:
390: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
392: Notes:
393: the initial solution and current solution are not display with a common axis scaling so generally the option -ts_monitor_draw_solution_initial
394: will look bad
396: Level: intermediate
398: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
399: @*/
400: PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
401: {
402: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
403: PetscDraw draw;
405: if (!step && ictx->showinitial) {
406: if (!ictx->initialsolution) {
407: VecDuplicate(u,&ictx->initialsolution);
408: }
409: VecCopy(u,ictx->initialsolution);
410: }
411: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return 0;
413: if (ictx->showinitial) {
414: PetscReal pause;
415: PetscViewerDrawGetPause(ictx->viewer,&pause);
416: PetscViewerDrawSetPause(ictx->viewer,0.0);
417: VecView(ictx->initialsolution,ictx->viewer);
418: PetscViewerDrawSetPause(ictx->viewer,pause);
419: PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
420: }
421: VecView(u,ictx->viewer);
422: if (ictx->showtimestepandtime) {
423: PetscReal xl,yl,xr,yr,h;
424: char time[32];
426: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
427: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
428: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
429: h = yl + .95*(yr - yl);
430: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
431: PetscDrawFlush(draw);
432: }
434: if (ictx->showinitial) {
435: PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
436: }
437: return 0;
438: }
440: /*@C
441: TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram
443: Collective on TS
445: Input Parameters:
446: + ts - the TS context
447: . step - current time-step
448: . ptime - current time
449: - dummy - either a viewer or NULL
451: Level: intermediate
453: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
454: @*/
455: PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
456: {
457: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
458: PetscDraw draw;
459: PetscDrawAxis axis;
460: PetscInt n;
461: PetscMPIInt size;
462: PetscReal U0,U1,xl,yl,xr,yr,h;
463: char time[32];
464: const PetscScalar *U;
465: PetscErrorCode ierr;
467: MPI_Comm_size(PetscObjectComm((PetscObject)ts),&size);
469: VecGetSize(u,&n);
472: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
473: PetscViewerDrawGetDrawAxis(ictx->viewer,0,&axis);
474: PetscDrawAxisGetLimits(axis,&xl,&xr,&yl,&yr);
475: if (!step) {
476: PetscDrawClear(draw);
477: PetscDrawAxisDraw(axis);
478: }
480: VecGetArrayRead(u,&U);
481: U0 = PetscRealPart(U[0]);
482: U1 = PetscRealPart(U[1]);
483: VecRestoreArrayRead(u,&U);
484: if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) return 0;
486: PetscDrawCollectiveBegin(draw);
487: PetscDrawPoint(draw,U0,U1,PETSC_DRAW_BLACK);
488: if (ictx->showtimestepandtime) {
489: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
490: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
491: h = yl + .95*(yr - yl);
492: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
493: }
494: PetscDrawCollectiveEnd(draw);
495: PetscDrawFlush(draw);
496: PetscDrawPause(draw);
497: PetscDrawSave(draw);
498: return 0;
499: }
501: /*@C
502: TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()
504: Collective on TS
506: Input Parameters:
507: . ctx - the monitor context
509: Level: intermediate
511: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
512: @*/
513: PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
514: {
515: PetscViewerDestroy(&(*ictx)->viewer);
516: VecDestroy(&(*ictx)->initialsolution);
517: PetscFree(*ictx);
518: return 0;
519: }
521: /*@C
522: TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx
524: Collective on TS
526: Input Parameter:
527: . ts - time-step context
529: Output Patameter:
530: . ctx - the monitor context
532: Options Database:
533: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
535: Level: intermediate
537: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
538: @*/
539: PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
540: {
541: PetscNew(ctx);
542: PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
543: PetscViewerSetFromOptions((*ctx)->viewer);
545: (*ctx)->howoften = howoften;
546: (*ctx)->showinitial = PETSC_FALSE;
547: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);
549: (*ctx)->showtimestepandtime = PETSC_FALSE;
550: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
551: return 0;
552: }
554: /*@C
555: TSMonitorDrawSolutionFunction - Monitors progress of the TS solvers by calling
556: VecView() for the solution provided by TSSetSolutionFunction() at each timestep
558: Collective on TS
560: Input Parameters:
561: + ts - the TS context
562: . step - current time-step
563: . ptime - current time
564: - dummy - either a viewer or NULL
566: Options Database:
567: . -ts_monitor_draw_solution_function - Monitor error graphically, requires user to have provided TSSetSolutionFunction()
569: Level: intermediate
571: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
572: @*/
573: PetscErrorCode TSMonitorDrawSolutionFunction(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
574: {
575: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
576: PetscViewer viewer = ctx->viewer;
577: Vec work;
579: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return 0;
580: VecDuplicate(u,&work);
581: TSComputeSolutionFunction(ts,ptime,work);
582: VecView(work,viewer);
583: VecDestroy(&work);
584: return 0;
585: }
587: /*@C
588: TSMonitorDrawError - Monitors progress of the TS solvers by calling
589: VecView() for the error at each timestep
591: Collective on TS
593: Input Parameters:
594: + ts - the TS context
595: . step - current time-step
596: . ptime - current time
597: - dummy - either a viewer or NULL
599: Options Database:
600: . -ts_monitor_draw_error - Monitor error graphically, requires user to have provided TSSetSolutionFunction()
602: Level: intermediate
604: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
605: @*/
606: PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
607: {
608: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
609: PetscViewer viewer = ctx->viewer;
610: Vec work;
612: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return 0;
613: VecDuplicate(u,&work);
614: TSComputeSolutionFunction(ts,ptime,work);
615: VecAXPY(work,-1.0,u);
616: VecView(work,viewer);
617: VecDestroy(&work);
618: return 0;
619: }
621: /*@C
622: TSMonitorSolution - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file or a PetscDraw object
624: Collective on TS
626: Input Parameters:
627: + ts - the TS context
628: . step - current time-step
629: . ptime - current time
630: . u - current state
631: - vf - viewer and its format
633: Level: intermediate
635: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
636: @*/
637: PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
638: {
639: PetscViewerPushFormat(vf->viewer,vf->format);
640: VecView(u,vf->viewer);
641: PetscViewerPopFormat(vf->viewer);
642: return 0;
643: }
645: /*@C
646: TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.
648: Collective on TS
650: Input Parameters:
651: + ts - the TS context
652: . step - current time-step
653: . ptime - current time
654: . u - current state
655: - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
657: Level: intermediate
659: Notes:
660: The VTK format does not allow writing multiple time steps in the same file, therefore a different file will be written for each time step.
661: These are named according to the file name template.
663: This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().
665: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
666: @*/
667: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
668: {
669: char filename[PETSC_MAX_PATH_LEN];
670: PetscViewer viewer;
672: if (step < 0) return 0; /* -1 indicates interpolated solution */
673: PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
674: PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
675: VecView(u,viewer);
676: PetscViewerDestroy(&viewer);
677: return 0;
678: }
680: /*@C
681: TSMonitorSolutionVTKDestroy - Destroy context for monitoring
683: Collective on TS
685: Input Parameters:
686: . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
688: Level: intermediate
690: Note:
691: This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().
693: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
694: @*/
695: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
696: {
697: PetscFree(*(char**)filenametemplate);
698: return 0;
699: }
701: /*@C
702: TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
703: in a time based line graph
705: Collective on TS
707: Input Parameters:
708: + ts - the TS context
709: . step - current time-step
710: . ptime - current time
711: . u - current solution
712: - dctx - the TSMonitorLGCtx object that contains all the options for the monitoring, this is created with TSMonitorLGCtxCreate()
714: Options Database:
715: . -ts_monitor_lg_solution_variables - enable monitor of lg solution variables
717: Level: intermediate
719: Notes:
720: Each process in a parallel run displays its component solutions in a separate window
722: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
723: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
724: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
725: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
726: @*/
727: PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
728: {
729: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dctx;
730: const PetscScalar *yy;
731: Vec v;
733: if (step < 0) return 0; /* -1 indicates interpolated solution */
734: if (!step) {
735: PetscDrawAxis axis;
736: PetscInt dim;
737: PetscDrawLGGetAxis(ctx->lg,&axis);
738: PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
739: if (!ctx->names) {
740: PetscBool flg;
741: /* user provides names of variables to plot but no names has been set so assume names are integer values */
742: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",&flg);
743: if (flg) {
744: PetscInt i,n;
745: char **names;
746: VecGetSize(u,&n);
747: PetscMalloc1(n+1,&names);
748: for (i=0; i<n; i++) {
749: PetscMalloc1(5,&names[i]);
750: PetscSNPrintf(names[i],5,"%D",i);
751: }
752: names[n] = NULL;
753: ctx->names = names;
754: }
755: }
756: if (ctx->names && !ctx->displaynames) {
757: char **displaynames;
758: PetscBool flg;
759: VecGetLocalSize(u,&dim);
760: PetscCalloc1(dim+1,&displaynames);
761: PetscOptionsGetStringArray(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);
762: if (flg) {
763: TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);
764: }
765: PetscStrArrayDestroy(&displaynames);
766: }
767: if (ctx->displaynames) {
768: PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);
769: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);
770: } else if (ctx->names) {
771: VecGetLocalSize(u,&dim);
772: PetscDrawLGSetDimension(ctx->lg,dim);
773: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);
774: } else {
775: VecGetLocalSize(u,&dim);
776: PetscDrawLGSetDimension(ctx->lg,dim);
777: }
778: PetscDrawLGReset(ctx->lg);
779: }
781: if (!ctx->transform) v = u;
782: else (*ctx->transform)(ctx->transformctx,u,&v);
783: VecGetArrayRead(v,&yy);
784: if (ctx->displaynames) {
785: PetscInt i;
786: for (i=0; i<ctx->ndisplayvariables; i++)
787: ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
788: PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);
789: } else {
790: #if defined(PETSC_USE_COMPLEX)
791: PetscInt i,n;
792: PetscReal *yreal;
793: VecGetLocalSize(v,&n);
794: PetscMalloc1(n,&yreal);
795: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
796: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
797: PetscFree(yreal);
798: #else
799: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
800: #endif
801: }
802: VecRestoreArrayRead(v,&yy);
803: if (ctx->transform) VecDestroy(&v);
805: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
806: PetscDrawLGDraw(ctx->lg);
807: PetscDrawLGSave(ctx->lg);
808: }
809: return 0;
810: }
812: /*@C
813: TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
815: Collective on TS
817: Input Parameters:
818: + ts - the TS context
819: - names - the names of the components, final string must be NULL
821: Level: intermediate
823: Notes:
824: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
826: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames()
827: @*/
828: PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names)
829: {
830: PetscInt i;
832: for (i=0; i<ts->numbermonitors; i++) {
833: if (ts->monitor[i] == TSMonitorLGSolution) {
834: TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);
835: break;
836: }
837: }
838: return 0;
839: }
841: /*@C
842: TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
844: Collective on TS
846: Input Parameters:
847: + ts - the TS context
848: - names - the names of the components, final string must be NULL
850: Level: intermediate
852: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames()
853: @*/
854: PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names)
855: {
856: PetscStrArrayDestroy(&ctx->names);
857: PetscStrArrayallocpy(names,&ctx->names);
858: return 0;
859: }
861: /*@C
862: TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot
864: Collective on TS
866: Input Parameter:
867: . ts - the TS context
869: Output Parameter:
870: . names - the names of the components, final string must be NULL
872: Level: intermediate
874: Notes:
875: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
877: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
878: @*/
879: PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names)
880: {
881: PetscInt i;
883: *names = NULL;
884: for (i=0; i<ts->numbermonitors; i++) {
885: if (ts->monitor[i] == TSMonitorLGSolution) {
886: TSMonitorLGCtx ctx = (TSMonitorLGCtx) ts->monitorcontext[i];
887: *names = (const char *const *)ctx->names;
888: break;
889: }
890: }
891: return 0;
892: }
894: /*@C
895: TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor
897: Collective on TS
899: Input Parameters:
900: + ctx - the TSMonitorLG context
901: - displaynames - the names of the components, final string must be NULL
903: Level: intermediate
905: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
906: @*/
907: PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames)
908: {
909: PetscInt j = 0,k;
911: if (!ctx->names) return 0;
912: PetscStrArrayDestroy(&ctx->displaynames);
913: PetscStrArrayallocpy(displaynames,&ctx->displaynames);
914: while (displaynames[j]) j++;
915: ctx->ndisplayvariables = j;
916: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);
917: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);
918: j = 0;
919: while (displaynames[j]) {
920: k = 0;
921: while (ctx->names[k]) {
922: PetscBool flg;
923: PetscStrcmp(displaynames[j],ctx->names[k],&flg);
924: if (flg) {
925: ctx->displayvariables[j] = k;
926: break;
927: }
928: k++;
929: }
930: j++;
931: }
932: return 0;
933: }
935: /*@C
936: TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor
938: Collective on TS
940: Input Parameters:
941: + ts - the TS context
942: - displaynames - the names of the components, final string must be NULL
944: Notes:
945: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
947: Level: intermediate
949: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
950: @*/
951: PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames)
952: {
953: PetscInt i;
955: for (i=0; i<ts->numbermonitors; i++) {
956: if (ts->monitor[i] == TSMonitorLGSolution) {
957: TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);
958: break;
959: }
960: }
961: return 0;
962: }
964: /*@C
965: TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed
967: Collective on TS
969: Input Parameters:
970: + ts - the TS context
971: . transform - the transform function
972: . destroy - function to destroy the optional context
973: - ctx - optional context used by transform function
975: Notes:
976: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
978: Level: intermediate
980: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform()
981: @*/
982: PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
983: {
984: PetscInt i;
986: for (i=0; i<ts->numbermonitors; i++) {
987: if (ts->monitor[i] == TSMonitorLGSolution) {
988: TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);
989: }
990: }
991: return 0;
992: }
994: /*@C
995: TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed
997: Collective on TSLGCtx
999: Input Parameters:
1000: + ts - the TS context
1001: . transform - the transform function
1002: . destroy - function to destroy the optional context
1003: - ctx - optional context used by transform function
1005: Level: intermediate
1007: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform()
1008: @*/
1009: PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
1010: {
1011: ctx->transform = transform;
1012: ctx->transformdestroy = destroy;
1013: ctx->transformctx = tctx;
1014: return 0;
1015: }
1017: /*@C
1018: TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the error
1019: in a time based line graph
1021: Collective on TS
1023: Input Parameters:
1024: + ts - the TS context
1025: . step - current time-step
1026: . ptime - current time
1027: . u - current solution
1028: - dctx - TSMonitorLGCtx object created with TSMonitorLGCtxCreate()
1030: Level: intermediate
1032: Notes:
1033: Each process in a parallel run displays its component errors in a separate window
1035: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
1037: Options Database Keys:
1038: . -ts_monitor_lg_error - create a graphical monitor of error history
1040: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
1041: @*/
1042: PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
1043: {
1044: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy;
1045: const PetscScalar *yy;
1046: Vec y;
1048: if (!step) {
1049: PetscDrawAxis axis;
1050: PetscInt dim;
1051: PetscDrawLGGetAxis(ctx->lg,&axis);
1052: PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Error");
1053: VecGetLocalSize(u,&dim);
1054: PetscDrawLGSetDimension(ctx->lg,dim);
1055: PetscDrawLGReset(ctx->lg);
1056: }
1057: VecDuplicate(u,&y);
1058: TSComputeSolutionFunction(ts,ptime,y);
1059: VecAXPY(y,-1.0,u);
1060: VecGetArrayRead(y,&yy);
1061: #if defined(PETSC_USE_COMPLEX)
1062: {
1063: PetscReal *yreal;
1064: PetscInt i,n;
1065: VecGetLocalSize(y,&n);
1066: PetscMalloc1(n,&yreal);
1067: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
1068: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
1069: PetscFree(yreal);
1070: }
1071: #else
1072: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
1073: #endif
1074: VecRestoreArrayRead(y,&yy);
1075: VecDestroy(&y);
1076: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
1077: PetscDrawLGDraw(ctx->lg);
1078: PetscDrawLGSave(ctx->lg);
1079: }
1080: return 0;
1081: }
1083: /*@C
1084: TSMonitorSPSwarmSolution - Graphically displays phase plots of DMSwarm particles on a scatter plot
1086: Input Parameters:
1087: + ts - the TS context
1088: . step - current time-step
1089: . ptime - current time
1090: . u - current solution
1091: - dctx - the TSMonitorSPCtx object that contains all the options for the monitoring, this is created with TSMonitorSPCtxCreate()
1093: Options Database:
1094: + -ts_monitor_sp_swarm <n> - Monitor the solution every n steps, or -1 for plotting only the final solution
1095: . -ts_monitor_sp_swarm_retain <n> - Retain n old points so we can see the history, or -1 for all points
1096: - -ts_monitor_sp_swarm_phase <bool> - Plot in phase space, as opposed to coordinate space
1098: Level: intermediate
1100: .seealso: TSMonitoSet()
1101: @*/
1102: PetscErrorCode TSMonitorSPSwarmSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dctx)
1103: {
1104: TSMonitorSPCtx ctx = (TSMonitorSPCtx) dctx;
1105: DM dm, cdm;
1106: const PetscScalar *yy;
1107: PetscReal *y, *x;
1108: PetscInt Np, p, dim = 2;
1110: if (step < 0) return 0; /* -1 indicates interpolated solution */
1111: if (!step) {
1112: PetscDrawAxis axis;
1113: PetscReal dmboxlower[2], dmboxupper[2];
1114: TSGetDM(ts, &dm);
1115: DMGetDimension(dm, &dim);
1117: DMSwarmGetCellDM(dm, &cdm);
1118: DMGetBoundingBox(cdm, dmboxlower, dmboxupper);
1119: VecGetLocalSize(u, &Np);
1120: Np /= dim*2;
1121: PetscDrawSPGetAxis(ctx->sp,&axis);
1122: if (ctx->phase) {
1123: PetscDrawAxisSetLabels(axis,"Particles","X","V");
1124: PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], -5, 5);
1125: } else {
1126: PetscDrawAxisSetLabels(axis,"Particles","X","Y");
1127: PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], dmboxlower[1], dmboxupper[1]);
1128: }
1129: PetscDrawAxisSetHoldLimits(axis, PETSC_TRUE);
1130: PetscDrawSPSetDimension(ctx->sp, Np);
1131: PetscDrawSPReset(ctx->sp);
1132: }
1133: VecGetLocalSize(u, &Np);
1134: Np /= dim*2;
1135: VecGetArrayRead(u,&yy);
1136: PetscMalloc2(Np, &x, Np, &y);
1137: /* get points from solution vector */
1138: for (p = 0; p < Np; ++p) {
1139: if (ctx->phase) {
1140: x[p] = PetscRealPart(yy[p*dim*2]);
1141: y[p] = PetscRealPart(yy[p*dim*2 + dim]);
1142: } else {
1143: x[p] = PetscRealPart(yy[p*dim*2]);
1144: y[p] = PetscRealPart(yy[p*dim*2 + 1]);
1145: }
1146: }
1147: VecRestoreArrayRead(u,&yy);
1148: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
1149: PetscDraw draw;
1150: PetscDrawSPGetDraw(ctx->sp, &draw);
1151: if ((ctx->retain == 0) || (ctx->retain > 0 && !(step % ctx->retain))) {
1152: PetscDrawClear(draw);
1153: }
1154: PetscDrawFlush(draw);
1155: PetscDrawSPReset(ctx->sp);
1156: PetscDrawSPAddPoint(ctx->sp, x, y);
1157: PetscDrawSPDraw(ctx->sp, PETSC_FALSE);
1158: PetscDrawSPSave(ctx->sp);
1159: }
1160: PetscFree2(x, y);
1161: return 0;
1162: }
1164: /*@C
1165: TSMonitorError - Monitors progress of the TS solvers by printing the 2 norm of the error at each timestep
1167: Collective on TS
1169: Input Parameters:
1170: + ts - the TS context
1171: . step - current time-step
1172: . ptime - current time
1173: . u - current solution
1174: - dctx - unused context
1176: Level: intermediate
1178: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
1180: Options Database Keys:
1181: . -ts_monitor_error - create a graphical monitor of error history
1183: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
1184: @*/
1185: PetscErrorCode TSMonitorError(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
1186: {
1187: DM dm;
1188: PetscDS ds = NULL;
1189: PetscInt Nf = -1, f;
1190: PetscBool flg;
1192: TSGetDM(ts, &dm);
1193: if (dm) DMGetDS(dm, &ds);
1194: if (ds) PetscDSGetNumFields(ds, &Nf);
1195: if (Nf <= 0) {
1196: Vec y;
1197: PetscReal nrm;
1199: VecDuplicate(u,&y);
1200: TSComputeSolutionFunction(ts,ptime,y);
1201: VecAXPY(y,-1.0,u);
1202: PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERASCII,&flg);
1203: if (flg) {
1204: VecNorm(y,NORM_2,&nrm);
1205: PetscViewerASCIIPrintf(vf->viewer,"2-norm of error %g\n",(double)nrm);
1206: }
1207: PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERDRAW,&flg);
1208: if (flg) {
1209: VecView(y,vf->viewer);
1210: }
1211: VecDestroy(&y);
1212: } else {
1213: PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
1214: void **ctxs;
1215: Vec v;
1216: PetscReal ferrors[1];
1218: PetscMalloc2(Nf, &exactFuncs, Nf, &ctxs);
1219: for (f = 0; f < Nf; ++f) PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]);
1220: DMComputeL2FieldDiff(dm, ptime, exactFuncs, ctxs, u, ferrors);
1221: PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [", (int) step, (double) ptime);
1222: for (f = 0; f < Nf; ++f) {
1223: if (f > 0) PetscPrintf(PETSC_COMM_WORLD, ", ");
1224: PetscPrintf(PETSC_COMM_WORLD, "%2.3g", (double) ferrors[f]);
1225: }
1226: PetscPrintf(PETSC_COMM_WORLD, "]\n");
1228: VecViewFromOptions(u, NULL, "-sol_vec_view");
1230: PetscOptionsHasName(NULL, NULL, "-exact_vec_view", &flg);
1231: if (flg) {
1232: DMGetGlobalVector(dm, &v);
1233: DMProjectFunction(dm, ptime, exactFuncs, ctxs, INSERT_ALL_VALUES, v);
1234: PetscObjectSetName((PetscObject) v, "Exact Solution");
1235: VecViewFromOptions(v, NULL, "-exact_vec_view");
1236: DMRestoreGlobalVector(dm, &v);
1237: }
1238: PetscFree2(exactFuncs, ctxs);
1239: }
1240: return 0;
1241: }
1243: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
1244: {
1245: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
1246: PetscReal x = ptime,y;
1247: PetscInt its;
1249: if (n < 0) return 0; /* -1 indicates interpolated solution */
1250: if (!n) {
1251: PetscDrawAxis axis;
1252: PetscDrawLGGetAxis(ctx->lg,&axis);
1253: PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
1254: PetscDrawLGReset(ctx->lg);
1255: ctx->snes_its = 0;
1256: }
1257: TSGetSNESIterations(ts,&its);
1258: y = its - ctx->snes_its;
1259: PetscDrawLGAddPoint(ctx->lg,&x,&y);
1260: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
1261: PetscDrawLGDraw(ctx->lg);
1262: PetscDrawLGSave(ctx->lg);
1263: }
1264: ctx->snes_its = its;
1265: return 0;
1266: }
1268: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
1269: {
1270: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
1271: PetscReal x = ptime,y;
1272: PetscInt its;
1274: if (n < 0) return 0; /* -1 indicates interpolated solution */
1275: if (!n) {
1276: PetscDrawAxis axis;
1277: PetscDrawLGGetAxis(ctx->lg,&axis);
1278: PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
1279: PetscDrawLGReset(ctx->lg);
1280: ctx->ksp_its = 0;
1281: }
1282: TSGetKSPIterations(ts,&its);
1283: y = its - ctx->ksp_its;
1284: PetscDrawLGAddPoint(ctx->lg,&x,&y);
1285: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
1286: PetscDrawLGDraw(ctx->lg);
1287: PetscDrawLGSave(ctx->lg);
1288: }
1289: ctx->ksp_its = its;
1290: return 0;
1291: }
1293: /*@C
1294: TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()
1296: Collective on TS
1298: Input Parameters:
1299: . ts - the ODE solver object
1301: Output Parameter:
1302: . ctx - the context
1304: Level: intermediate
1306: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
1308: @*/
1309: PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
1310: {
1311: PetscNew(ctx);
1312: return 0;
1313: }
1315: /*@C
1316: TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution
1318: Collective on TS
1320: Input Parameters:
1321: + ts - the TS context
1322: . step - current time-step
1323: . ptime - current time
1324: . u - current solution
1325: - dctx - the envelope context
1327: Options Database:
1328: . -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
1330: Level: intermediate
1332: Notes:
1333: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope
1335: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxCreate()
1336: @*/
1337: PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
1338: {
1339: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;
1341: if (!ctx->max) {
1342: VecDuplicate(u,&ctx->max);
1343: VecDuplicate(u,&ctx->min);
1344: VecCopy(u,ctx->max);
1345: VecCopy(u,ctx->min);
1346: } else {
1347: VecPointwiseMax(ctx->max,u,ctx->max);
1348: VecPointwiseMin(ctx->min,u,ctx->min);
1349: }
1350: return 0;
1351: }
1353: /*@C
1354: TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution
1356: Collective on TS
1358: Input Parameter:
1359: . ts - the TS context
1361: Output Parameters:
1362: + max - the maximum values
1363: - min - the minimum values
1365: Notes:
1366: If the TS does not have a TSMonitorEnvelopeCtx associated with it then this function is ignored
1368: Level: intermediate
1370: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
1371: @*/
1372: PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
1373: {
1374: PetscInt i;
1376: if (max) *max = NULL;
1377: if (min) *min = NULL;
1378: for (i=0; i<ts->numbermonitors; i++) {
1379: if (ts->monitor[i] == TSMonitorEnvelope) {
1380: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i];
1381: if (max) *max = ctx->max;
1382: if (min) *min = ctx->min;
1383: break;
1384: }
1385: }
1386: return 0;
1387: }
1389: /*@C
1390: TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate().
1392: Collective on TSMonitorEnvelopeCtx
1394: Input Parameter:
1395: . ctx - the monitor context
1397: Level: intermediate
1399: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep()
1400: @*/
1401: PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
1402: {
1403: VecDestroy(&(*ctx)->min);
1404: VecDestroy(&(*ctx)->max);
1405: PetscFree(*ctx);
1406: return 0;
1407: }
1409: /*@C
1410: TSDMSwarmMonitorMoments - Monitors the first three moments of a DMSarm being evolved by the TS
1412: Not collective
1414: Input Parameters:
1415: + ts - the TS context
1416: . step - current timestep
1417: . t - current time
1418: . u - current solution
1419: - ctx - not used
1421: Options Database:
1422: . -ts_dmswarm_monitor_moments - Monitor moments of particle distribution
1424: Level: intermediate
1426: Notes:
1427: This requires a DMSwarm be attached to the TS.
1429: .seealso: TSMonitorSet(), TSMonitorDefault(), DMSWARM
1430: @*/
1431: PetscErrorCode TSDMSwarmMonitorMoments(TS ts, PetscInt step, PetscReal t, Vec U, PetscViewerAndFormat *vf)
1432: {
1433: DM sw;
1434: const PetscScalar *u;
1435: PetscReal m = 1.0, totE = 0., totMom[3] = {0., 0., 0.};
1436: PetscInt dim, d, Np, p;
1437: MPI_Comm comm;
1440: TSGetDM(ts, &sw);
1441: if (!sw || step%ts->monitorFrequency != 0) return 0;
1442: PetscObjectGetComm((PetscObject) ts, &comm);
1443: DMGetDimension(sw, &dim);
1444: VecGetLocalSize(U, &Np);
1445: Np /= dim;
1446: VecGetArrayRead(U, &u);
1447: for (p = 0; p < Np; ++p) {
1448: for (d = 0; d < dim; ++d) {
1449: totE += PetscRealPart(u[p*dim+d]*u[p*dim+d]);
1450: totMom[d] += PetscRealPart(u[p*dim+d]);
1451: }
1452: }
1453: VecRestoreArrayRead(U, &u);
1454: for (d = 0; d < dim; ++d) totMom[d] *= m;
1455: totE *= 0.5*m;
1456: PetscPrintf(comm, "Step %4D Total Energy: %10.8lf", step, (double) totE);
1457: for (d = 0; d < dim; ++d) PetscPrintf(comm, " Total Momentum %c: %10.8lf", 'x'+d, (double) totMom[d]);
1458: PetscPrintf(comm, "\n");
1459: return 0;
1460: }