xref: /petsc/src/ts/interface/tsmon.c (revision bfe80ac4a46d58cb7760074b25f5e81b2f541d8a)
1 #include <petsc/private/tsimpl.h> /*I "petscts.h"  I*/
2 #include <petscdm.h>
3 #include <petscds.h>
4 #include <petscdmswarm.h>
5 #include <petscdraw.h>
6 
7 /*@C
8   TSMonitor - Runs all user-provided monitor routines set using `TSMonitorSet()`
9 
10   Collective
11 
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
17 
18   Level: developer
19 
20   Notes:
21   `TSMonitor()` is typically used automatically within the time stepping implementations.
22   Users would almost never call this routine directly.
23 
24   A step of -1 indicates that the monitor is being called on a solution obtained by interpolating from computed solutions
25 
26 .seealso: `TS`, `TSMonitorSet()`, `TSMonitorSetFromOptions()`
27 @*/
28 PetscErrorCode TSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec u)
29 {
30   DM       dm;
31   PetscInt i, n = ts->numbermonitors;
32 
33   PetscFunctionBegin;
34   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
35   PetscValidHeaderSpecific(u, VEC_CLASSID, 4);
36 
37   PetscCall(TSGetDM(ts, &dm));
38   PetscCall(DMSetOutputSequenceNumber(dm, step, ptime));
39 
40   PetscCall(VecLockReadPush(u));
41   for (i = 0; i < n; i++) PetscCall((*ts->monitor[i])(ts, step, ptime, u, ts->monitorcontext[i]));
42   PetscCall(VecLockReadPop(u));
43   PetscFunctionReturn(PETSC_SUCCESS);
44 }
45 
46 /*@C
47   TSMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
48 
49   Collective
50 
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, this must use a `PetscViewerFormat` as its context
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
58 
59   Level: developer
60 
61 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `PetscOptionsCreateViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
62           `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
63           `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`,
64           `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
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;
74 
75   PetscFunctionBegin;
76   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, name, &viewer, &format, &flg));
77   if (flg) {
78     PetscViewerAndFormat *vf;
79     char                  interval_key[1024];
80 
81     PetscCall(PetscSNPrintf(interval_key, sizeof interval_key, "%s_interval", name));
82     PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf));
83     vf->view_interval = 1;
84     PetscCall(PetscOptionsGetInt(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, interval_key, &vf->view_interval, NULL));
85     PetscCall(PetscViewerDestroy(&viewer));
86     if (monitorsetup) PetscCall((*monitorsetup)(ts, vf));
87     PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))monitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
88   }
89   PetscFunctionReturn(PETSC_SUCCESS);
90 }
91 
92 /*@C
93   TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
94   timestep to display the iteration's  progress.
95 
96   Logically Collective
97 
98   Input Parameters:
99 + ts       - the `TS` context obtained from `TSCreate()`
100 . monitor  - monitoring routine
101 . mctx     - [optional] user-defined context for private data for the monitor routine (use `NULL` if no context is desired)
102 - mdestroy - [optional] routine that frees monitor context (may be `NULL`), see `PetscCtxDestroyFn` for the calling sequence
103 
104   Calling sequence of `monitor`:
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 - ctx   - [optional] monitoring context
110 
111   Level: intermediate
112 
113   Note:
114   This routine adds an additional monitor to the list of monitors that already has been loaded.
115 
116   Fortran Notes:
117   Only a single monitor function can be set for each `TS` object
118 
119 .seealso: [](ch_ts), `TSMonitorDefault()`, `TSMonitorCancel()`, `TSDMSwarmMonitorMoments()`, `TSMonitorExtreme()`, `TSMonitorDrawSolution()`,
120           `TSMonitorDrawSolutionPhase()`, `TSMonitorDrawSolutionFunction()`, `TSMonitorDrawError()`, `TSMonitorSolution()`, `TSMonitorSolutionVTK()`,
121           `TSMonitorLGSolution()`, `TSMonitorLGError()`, `TSMonitorSPSwarmSolution()`, `TSMonitorError()`, `TSMonitorEnvelope()`,  `PetscCtxDestroyFn`
122 @*/
123 PetscErrorCode TSMonitorSet(TS ts, PetscErrorCode (*monitor)(TS ts, PetscInt steps, PetscReal time, Vec u, void *ctx), void *mctx, PetscCtxDestroyFn *mdestroy)
124 {
125   PetscInt  i;
126   PetscBool identical;
127 
128   PetscFunctionBegin;
129   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
130   for (i = 0; i < ts->numbermonitors; i++) {
131     PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))monitor, mctx, mdestroy, (PetscErrorCode (*)(void))ts->monitor[i], ts->monitorcontext[i], ts->monitordestroy[i], &identical));
132     if (identical) PetscFunctionReturn(PETSC_SUCCESS);
133   }
134   PetscCheck(ts->numbermonitors < MAXTSMONITORS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many monitors set");
135   ts->monitor[ts->numbermonitors]          = monitor;
136   ts->monitordestroy[ts->numbermonitors]   = mdestroy;
137   ts->monitorcontext[ts->numbermonitors++] = mctx;
138   PetscFunctionReturn(PETSC_SUCCESS);
139 }
140 
141 /*@C
142   TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
143 
144   Logically Collective
145 
146   Input Parameter:
147 . ts - the `TS` context obtained from `TSCreate()`
148 
149   Level: intermediate
150 
151   Note:
152   There is no way to remove a single, specific monitor.
153 
154 .seealso: [](ch_ts), `TS`, `TSMonitorDefault()`, `TSMonitorSet()`
155 @*/
156 PetscErrorCode TSMonitorCancel(TS ts)
157 {
158   PetscInt i;
159 
160   PetscFunctionBegin;
161   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
162   for (i = 0; i < ts->numbermonitors; i++) {
163     if (ts->monitordestroy[i]) PetscCall((*ts->monitordestroy[i])(&ts->monitorcontext[i]));
164   }
165   ts->numbermonitors = 0;
166   PetscFunctionReturn(PETSC_SUCCESS);
167 }
168 
169 /*@C
170   TSMonitorDefault - The default monitor, prints the timestep and time for each step
171 
172   Input Parameters:
173 + ts    - the `TS` context
174 . step  - 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)
175 . ptime - current time
176 . v     - current iterate
177 - vf    - the viewer and format
178 
179   Options Database Key:
180 . -ts_monitor - monitors the time integration
181 
182   Level: intermediate
183 
184   Notes:
185   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
186   to be used during the `TS` integration.
187 
188 .seealso: [](ch_ts), `TSMonitorSet()`, `TSDMSwarmMonitorMoments()`, `TSMonitorExtreme()`, `TSMonitorDrawSolution()`,
189           `TSMonitorDrawSolutionPhase()`, `TSMonitorDrawSolutionFunction()`, `TSMonitorDrawError()`, `TSMonitorSolution()`, `TSMonitorSolutionVTK()`,
190           `TSMonitorLGSolution()`, `TSMonitorLGError()`, `TSMonitorSPSwarmSolution()`, `TSMonitorError()`, `TSMonitorEnvelope()`
191 @*/
192 PetscErrorCode TSMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscViewerAndFormat *vf)
193 {
194   PetscViewer viewer = vf->viewer;
195   PetscBool   iascii, ibinary;
196 
197   PetscFunctionBegin;
198   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 5);
199   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
200   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
201   PetscCall(PetscViewerPushFormat(viewer, vf->format));
202   if (iascii) {
203     PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
204     if (step == -1) { /* this indicates it is an interpolated solution */
205       PetscCall(PetscViewerASCIIPrintf(viewer, "Interpolated solution at time %g between steps %" PetscInt_FMT " and %" PetscInt_FMT "\n", (double)ptime, ts->steps - 1, ts->steps));
206     } else {
207       PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " TS dt %g time %g%s", step, (double)ts->time_step, (double)ptime, ts->steprollback ? " (r)\n" : "\n"));
208     }
209     PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
210   } else if (ibinary) {
211     PetscMPIInt rank;
212     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
213     if (rank == 0) {
214       PetscBool skipHeader;
215       PetscInt  classid = REAL_FILE_CLASSID;
216 
217       PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
218       if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
219       PetscCall(PetscRealView(1, &ptime, viewer));
220     } else {
221       PetscCall(PetscRealView(0, &ptime, viewer));
222     }
223   }
224   PetscCall(PetscViewerPopFormat(viewer));
225   PetscFunctionReturn(PETSC_SUCCESS);
226 }
227 
228 /*@C
229   TSMonitorExtreme - Prints the extreme values of the solution at each timestep
230 
231   Input Parameters:
232 + ts    - the `TS` context
233 . step  - 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)
234 . ptime - current time
235 . v     - current iterate
236 - vf    - the viewer and format
237 
238   Level: intermediate
239 
240   Note:
241   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
242   to be used during the `TS` integration.
243 
244 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`
245 @*/
246 PetscErrorCode TSMonitorExtreme(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscViewerAndFormat *vf)
247 {
248   PetscViewer viewer = vf->viewer;
249   PetscBool   iascii;
250   PetscReal   max, min;
251 
252   PetscFunctionBegin;
253   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 5);
254   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
255   PetscCall(PetscViewerPushFormat(viewer, vf->format));
256   if (iascii) {
257     PetscCall(VecMax(v, NULL, &max));
258     PetscCall(VecMin(v, NULL, &min));
259     PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
260     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " 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));
261     PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
262   }
263   PetscCall(PetscViewerPopFormat(viewer));
264   PetscFunctionReturn(PETSC_SUCCESS);
265 }
266 
267 /*@C
268   TSMonitorLGCtxCreate - Creates a `TSMonitorLGCtx` context for use with
269   `TS` to monitor the solution process graphically in various ways
270 
271   Collective
272 
273   Input Parameters:
274 + comm     - the MPI communicator to use
275 . host     - the X display to open, or `NULL` for the local machine
276 . label    - the title to put in the title bar
277 . x        - the x screen coordinates of the upper left coordinate of the window
278 . y        - the y screen coordinates of the upper left coordinate of the window
279 . m        - the screen width in pixels
280 . n        - the screen height in pixels
281 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
282 
283   Output Parameter:
284 . ctx - the context
285 
286   Options Database Keys:
287 + -ts_monitor_lg_timestep        - automatically sets line graph monitor
288 . -ts_monitor_lg_timestep_log    - automatically sets line graph monitor
289 . -ts_monitor_lg_solution        - monitor the solution (or certain values of the solution by calling `TSMonitorLGSetDisplayVariables()` or `TSMonitorLGCtxSetDisplayVariables()`)
290 . -ts_monitor_lg_error           - monitor the error
291 . -ts_monitor_lg_ksp_iterations  - monitor the number of `KSP` iterations needed for each timestep
292 . -ts_monitor_lg_snes_iterations - monitor the number of `SNES` iterations needed for each timestep
293 - -lg_use_markers <true,false>   - mark the data points (at each time step) on the plot; default is true
294 
295   Level: intermediate
296 
297   Notes:
298   Pass the context and `TSMonitorLGCtxDestroy()` to `TSMonitorSet()` to have the context destroyed when no longer needed.
299 
300   One can provide a function that transforms the solution before plotting it with `TSMonitorLGCtxSetTransform()` or `TSMonitorLGSetTransform()`
301 
302   Many of the functions that control the monitoring have two forms\: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a `TS` object as the
303   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
304   as the first argument.
305 
306   One can control the names displayed for each solution or error variable with `TSMonitorLGCtxSetVariableNames()` or `TSMonitorLGSetVariableNames()`
307 
308 .seealso: [](ch_ts), `TSMonitorLGTimeStep()`, `TSMonitorSet()`, `TSMonitorLGSolution()`, `TSMonitorLGError()`, `TSMonitorDefault()`, `VecView()`,
309           `TSMonitorLGCtxSetVariableNames()`, `TSMonitorLGCtxGetVariableNames()`,
310           `TSMonitorLGSetVariableNames()`, `TSMonitorLGGetVariableNames()`, `TSMonitorLGSetDisplayVariables()`, `TSMonitorLGCtxSetDisplayVariables()`,
311           `TSMonitorLGCtxSetTransform()`, `TSMonitorLGSetTransform()`, `TSMonitorLGSNESIterations()`, `TSMonitorLGKSPIterations()`,
312           `TSMonitorEnvelopeCtxCreate()`, `TSMonitorEnvelopeGetBounds()`, `TSMonitorEnvelopeCtxDestroy()`, `TSMonitorEnvelop()`
313 @*/
314 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm, const char host[], const char label[], int x, int y, int m, int n, PetscInt howoften, TSMonitorLGCtx *ctx)
315 {
316   PetscDraw draw;
317 
318   PetscFunctionBegin;
319   PetscCall(PetscNew(ctx));
320   PetscCall(PetscDrawCreate(comm, host, label, x, y, m, n, &draw));
321   PetscCall(PetscDrawSetFromOptions(draw));
322   PetscCall(PetscDrawLGCreate(draw, 1, &(*ctx)->lg));
323   PetscCall(PetscDrawLGSetFromOptions((*ctx)->lg));
324   PetscCall(PetscDrawDestroy(&draw));
325   (*ctx)->howoften = howoften;
326   PetscFunctionReturn(PETSC_SUCCESS);
327 }
328 
329 /*@C
330   TSMonitorLGTimeStep - Monitors a `TS` by printing the time-steps
331 
332   Collective
333 
334   Input Parameters:
335 + ts     - the time integrator
336 . step   - the current time step
337 . ptime  - the current time
338 . v      - the current state
339 - monctx - the monitor context obtained with `TSMonitorLGCtxCreate()`
340 
341   Level: advanced
342 
343   Note:
344   This is not called directly by users, rather one calls `TSMonitorSet()` along the `ctx` created by `TSMonitorLGCtxCreate()`
345   and `TSMonitorLGCtxDestroy()`
346 
347 .seealso: [](ch_ts), `TS`, `TSMonitorLGCtxCreate()`, `TSMonitorSet()`, `TSMonitorLGCtxDestroy()`
348 @*/
349 PetscErrorCode TSMonitorLGTimeStep(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx)
350 {
351   TSMonitorLGCtx ctx = (TSMonitorLGCtx)monctx;
352   PetscReal      x   = ptime, y;
353 
354   PetscFunctionBegin;
355   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates an interpolated solution */
356   if (!step) {
357     PetscDrawAxis axis;
358     const char   *ylabel = ctx->semilogy ? "Log Time Step" : "Time Step";
359     PetscCall(PetscDrawLGGetAxis(ctx->lg, &axis));
360     PetscCall(PetscDrawAxisSetLabels(axis, "Timestep as function of time", "Time", ylabel));
361     PetscCall(PetscDrawLGReset(ctx->lg));
362   }
363   PetscCall(TSGetTimeStep(ts, &y));
364   if (ctx->semilogy) y = PetscLog10Real(y);
365   PetscCall(PetscDrawLGAddPoint(ctx->lg, &x, &y));
366   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
367     PetscCall(PetscDrawLGDraw(ctx->lg));
368     PetscCall(PetscDrawLGSave(ctx->lg));
369   }
370   PetscFunctionReturn(PETSC_SUCCESS);
371 }
372 
373 /*@C
374   TSMonitorLGCtxDestroy - Destroys a line graph context that was created with `TSMonitorLGCtxCreate()`.
375 
376   Collective
377 
378   Input Parameter:
379 . ctx - the monitor context
380 
381   Level: intermediate
382 
383   Note:
384   Pass to `TSMonitorSet()` along with the context and `TSMonitorLGTimeStep()`
385 
386 .seealso: [](ch_ts), `TS`, `TSMonitorLGCtxCreate()`, `TSMonitorSet()`, `TSMonitorLGTimeStep()`
387 @*/
388 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
389 {
390   PetscFunctionBegin;
391   if ((*ctx)->transformdestroy) PetscCall(((*ctx)->transformdestroy)((void **)&(*ctx)->transformctx));
392   PetscCall(PetscDrawLGDestroy(&(*ctx)->lg));
393   PetscCall(PetscStrArrayDestroy(&(*ctx)->names));
394   PetscCall(PetscStrArrayDestroy(&(*ctx)->displaynames));
395   PetscCall(PetscFree((*ctx)->displayvariables));
396   PetscCall(PetscFree((*ctx)->displayvalues));
397   PetscCall(PetscFree(*ctx));
398   PetscFunctionReturn(PETSC_SUCCESS);
399 }
400 
401 /* Creates a TSMonitorSPCtx for use with DMSwarm particle visualizations */
402 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, PetscBool multispecies, TSMonitorSPCtx *ctx)
403 {
404   PetscDraw draw;
405 
406   PetscFunctionBegin;
407   PetscCall(PetscNew(ctx));
408   PetscCall(PetscDrawCreate(comm, host, label, x, y, m, n, &draw));
409   PetscCall(PetscDrawSetFromOptions(draw));
410   PetscCall(PetscDrawSPCreate(draw, 1, &(*ctx)->sp));
411   PetscCall(PetscDrawDestroy(&draw));
412   (*ctx)->howoften     = howoften;
413   (*ctx)->retain       = retain;
414   (*ctx)->phase        = phase;
415   (*ctx)->multispecies = multispecies;
416   PetscFunctionReturn(PETSC_SUCCESS);
417 }
418 
419 /* Destroys a TSMonitorSPCtx that was created with TSMonitorSPCtxCreate */
420 PetscErrorCode TSMonitorSPCtxDestroy(TSMonitorSPCtx *ctx)
421 {
422   PetscFunctionBegin;
423   PetscCall(PetscDrawSPDestroy(&(*ctx)->sp));
424   PetscCall(PetscFree(*ctx));
425   PetscFunctionReturn(PETSC_SUCCESS);
426 }
427 
428 /* Creates a TSMonitorHGCtx for use with DMSwarm particle visualizations */
429 PetscErrorCode TSMonitorHGCtxCreate(MPI_Comm comm, const char host[], const char label[], int x, int y, int m, int n, PetscInt howoften, PetscInt Ns, PetscInt Nb, PetscBool velocity, TSMonitorHGCtx *ctx)
430 {
431   PetscDraw draw;
432   int       Nsi, Nbi;
433 
434   PetscFunctionBegin;
435   PetscCall(PetscMPIIntCast(Ns, &Nsi));
436   PetscCall(PetscMPIIntCast(Nb, &Nbi));
437   PetscCall(PetscNew(ctx));
438   PetscCall(PetscMalloc1(Ns, &(*ctx)->hg));
439   for (int s = 0; s < Nsi; ++s) {
440     PetscCall(PetscDrawCreate(comm, host, label, x + s * m, y, m, n, &draw));
441     PetscCall(PetscDrawSetFromOptions(draw));
442     PetscCall(PetscDrawHGCreate(draw, Nbi, &(*ctx)->hg[s]));
443     PetscCall(PetscDrawHGCalcStats((*ctx)->hg[s], PETSC_TRUE));
444     PetscCall(PetscDrawDestroy(&draw));
445   }
446   (*ctx)->howoften = howoften;
447   (*ctx)->Ns       = Ns;
448   (*ctx)->velocity = velocity;
449   PetscFunctionReturn(PETSC_SUCCESS);
450 }
451 
452 /* Destroys a TSMonitorHGCtx that was created with TSMonitorHGCtxCreate */
453 PetscErrorCode TSMonitorHGCtxDestroy(TSMonitorHGCtx *ctx)
454 {
455   PetscInt s;
456 
457   PetscFunctionBegin;
458   for (s = 0; s < (*ctx)->Ns; ++s) PetscCall(PetscDrawHGDestroy(&(*ctx)->hg[s]));
459   PetscCall(PetscFree((*ctx)->hg));
460   PetscCall(PetscFree(*ctx));
461   PetscFunctionReturn(PETSC_SUCCESS);
462 }
463 
464 /*@C
465   TSMonitorDrawSolution - Monitors progress of the `TS` solvers by calling
466   `VecView()` for the solution at each timestep
467 
468   Collective
469 
470   Input Parameters:
471 + ts    - the `TS` context
472 . step  - current time-step
473 . ptime - current time
474 . u     - the solution at the current time
475 - dummy - either a viewer or `NULL`
476 
477   Options Database Keys:
478 + -ts_monitor_draw_solution         - draw the solution at each time-step
479 - -ts_monitor_draw_solution_initial - show initial solution as well as current solution
480 
481   Level: intermediate
482 
483   Notes:
484   The initial solution and current solution are not displayed with a common axis scaling so generally the option `-ts_monitor_draw_solution_initial`
485   will look bad
486 
487   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, as well as the context created with
488   `TSMonitorDrawCtxCreate()` and the function `TSMonitorDrawCtxDestroy()` to cause the monitor to be used during the `TS` integration.
489 
490 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorDrawCtxCreate()`, `TSMonitorDrawCtxDestroy()`
491 @*/
492 PetscErrorCode TSMonitorDrawSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dummy)
493 {
494   TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
495   PetscDraw        draw;
496 
497   PetscFunctionBegin;
498   if (!step && ictx->showinitial) {
499     if (!ictx->initialsolution) PetscCall(VecDuplicate(u, &ictx->initialsolution));
500     PetscCall(VecCopy(u, ictx->initialsolution));
501   }
502   if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(PETSC_SUCCESS);
503 
504   if (ictx->showinitial) {
505     PetscReal pause;
506     PetscCall(PetscViewerDrawGetPause(ictx->viewer, &pause));
507     PetscCall(PetscViewerDrawSetPause(ictx->viewer, 0.0));
508     PetscCall(VecView(ictx->initialsolution, ictx->viewer));
509     PetscCall(PetscViewerDrawSetPause(ictx->viewer, pause));
510     PetscCall(PetscViewerDrawSetHold(ictx->viewer, PETSC_TRUE));
511   }
512   PetscCall(VecView(u, ictx->viewer));
513   if (ictx->showtimestepandtime) {
514     PetscReal xl, yl, xr, yr, h;
515     char      time[32];
516 
517     PetscCall(PetscViewerDrawGetDraw(ictx->viewer, 0, &draw));
518     PetscCall(PetscSNPrintf(time, 32, "Timestep %" PetscInt_FMT " Time %g", step, (double)ptime));
519     PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
520     h = yl + .95 * (yr - yl);
521     PetscCall(PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time));
522     PetscCall(PetscDrawFlush(draw));
523   }
524 
525   if (ictx->showinitial) PetscCall(PetscViewerDrawSetHold(ictx->viewer, PETSC_FALSE));
526   PetscFunctionReturn(PETSC_SUCCESS);
527 }
528 
529 /*@C
530   TSMonitorDrawSolutionPhase - Monitors progress of the `TS` solvers by plotting the solution as a phase diagram
531 
532   Collective
533 
534   Input Parameters:
535 + ts    - the `TS` context
536 . step  - current time-step
537 . ptime - current time
538 . u     - the solution at the current time
539 - dummy - either a viewer or `NULL`
540 
541   Level: intermediate
542 
543   Notes:
544   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
545   to be used during the `TS` integration.
546 
547 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`
548 @*/
549 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dummy)
550 {
551   TSMonitorDrawCtx   ictx = (TSMonitorDrawCtx)dummy;
552   PetscDraw          draw;
553   PetscDrawAxis      axis;
554   PetscInt           n;
555   PetscMPIInt        size;
556   PetscReal          U0, U1, xl, yl, xr, yr, h;
557   char               time[32];
558   const PetscScalar *U;
559 
560   PetscFunctionBegin;
561   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ts), &size));
562   PetscCheck(size == 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Only allowed for sequential runs");
563   PetscCall(VecGetSize(u, &n));
564   PetscCheck(n == 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Only for ODEs with two unknowns");
565 
566   PetscCall(PetscViewerDrawGetDraw(ictx->viewer, 0, &draw));
567   PetscCall(PetscViewerDrawGetDrawAxis(ictx->viewer, 0, &axis));
568   PetscCall(PetscDrawAxisGetLimits(axis, &xl, &xr, &yl, &yr));
569   if (!step) {
570     PetscCall(PetscDrawClear(draw));
571     PetscCall(PetscDrawAxisDraw(axis));
572   }
573 
574   PetscCall(VecGetArrayRead(u, &U));
575   U0 = PetscRealPart(U[0]);
576   U1 = PetscRealPart(U[1]);
577   PetscCall(VecRestoreArrayRead(u, &U));
578   if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) PetscFunctionReturn(PETSC_SUCCESS);
579 
580   PetscDrawCollectiveBegin(draw);
581   PetscCall(PetscDrawPoint(draw, U0, U1, PETSC_DRAW_BLACK));
582   if (ictx->showtimestepandtime) {
583     PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
584     PetscCall(PetscSNPrintf(time, 32, "Timestep %" PetscInt_FMT " Time %g", step, (double)ptime));
585     h = yl + .95 * (yr - yl);
586     PetscCall(PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time));
587   }
588   PetscDrawCollectiveEnd(draw);
589   PetscCall(PetscDrawFlush(draw));
590   PetscCall(PetscDrawPause(draw));
591   PetscCall(PetscDrawSave(draw));
592   PetscFunctionReturn(PETSC_SUCCESS);
593 }
594 
595 /*@C
596   TSMonitorDrawCtxDestroy - Destroys the monitor context for `TSMonitorDrawSolution()`
597 
598   Collective
599 
600   Input Parameter:
601 . ictx - the monitor context
602 
603   Level: intermediate
604 
605 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorDrawSolution()`, `TSMonitorDrawError()`, `TSMonitorDrawCtx`
606 @*/
607 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
608 {
609   PetscFunctionBegin;
610   PetscCall(PetscViewerDestroy(&(*ictx)->viewer));
611   PetscCall(VecDestroy(&(*ictx)->initialsolution));
612   PetscCall(PetscFree(*ictx));
613   PetscFunctionReturn(PETSC_SUCCESS);
614 }
615 
616 /*@C
617   TSMonitorDrawCtxCreate - Creates the monitor context for `TSMonitorDrawCtx`
618 
619   Collective
620 
621   Input Parameters:
622 + comm     - the MPI communicator to use
623 . host     - the X display to open, or `NULL` for the local machine
624 . label    - the title to put in the title bar
625 . x        - the x screen coordinates of the upper left coordinate of the window
626 . y        - the y screen coordinates of the upper left coordinate of the window
627 . m        - the screen width in pixels
628 . n        - the screen height in pixels
629 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
630 
631   Output Parameter:
632 . ctx - the monitor context
633 
634   Options Database Keys:
635 + -ts_monitor_draw_solution         - draw the solution at each time-step
636 - -ts_monitor_draw_solution_initial - show initial solution as well as current solution
637 
638   Level: intermediate
639 
640   Note:
641   The context created by this function, `PetscMonitorDrawSolution()`, and `TSMonitorDrawCtxDestroy()` should be passed together to `TSMonitorSet()`.
642 
643 .seealso: [](ch_ts), `TS`, `TSMonitorDrawCtxDestroy()`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorDrawCtx`, `PetscMonitorDrawSolution()`
644 @*/
645 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm, const char host[], const char label[], int x, int y, int m, int n, PetscInt howoften, TSMonitorDrawCtx *ctx)
646 {
647   PetscFunctionBegin;
648   PetscCall(PetscNew(ctx));
649   PetscCall(PetscViewerDrawOpen(comm, host, label, x, y, m, n, &(*ctx)->viewer));
650   PetscCall(PetscViewerSetFromOptions((*ctx)->viewer));
651 
652   (*ctx)->howoften    = howoften;
653   (*ctx)->showinitial = PETSC_FALSE;
654   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_draw_solution_initial", &(*ctx)->showinitial, NULL));
655 
656   (*ctx)->showtimestepandtime = PETSC_FALSE;
657   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_draw_solution_show_time", &(*ctx)->showtimestepandtime, NULL));
658   PetscFunctionReturn(PETSC_SUCCESS);
659 }
660 
661 /*@C
662   TSMonitorDrawSolutionFunction - Monitors progress of the `TS` solvers by calling
663   `VecView()` for the solution provided by `TSSetSolutionFunction()` at each timestep
664 
665   Collective
666 
667   Input Parameters:
668 + ts    - the `TS` context
669 . step  - current time-step
670 . ptime - current time
671 . u     - solution at current time
672 - dummy - either a viewer or `NULL`
673 
674   Options Database Key:
675 . -ts_monitor_draw_solution_function - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
676 
677   Level: intermediate
678 
679   Note:
680   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
681   to be used during the `TS` integration.
682 
683 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSSetSolutionFunction()`
684 @*/
685 PetscErrorCode TSMonitorDrawSolutionFunction(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dummy)
686 {
687   TSMonitorDrawCtx ctx    = (TSMonitorDrawCtx)dummy;
688   PetscViewer      viewer = ctx->viewer;
689   Vec              work;
690 
691   PetscFunctionBegin;
692   if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(PETSC_SUCCESS);
693   PetscCall(VecDuplicate(u, &work));
694   PetscCall(TSComputeSolutionFunction(ts, ptime, work));
695   PetscCall(VecView(work, viewer));
696   PetscCall(VecDestroy(&work));
697   PetscFunctionReturn(PETSC_SUCCESS);
698 }
699 
700 /*@C
701   TSMonitorDrawError - Monitors progress of the `TS` solvers by calling
702   `VecView()` for the error at each timestep
703 
704   Collective
705 
706   Input Parameters:
707 + ts    - the `TS` context
708 . step  - current time-step
709 . ptime - current time
710 . u     - solution at current time
711 - dummy - either a viewer or `NULL`
712 
713   Options Database Key:
714 . -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
715 
716   Level: intermediate
717 
718   Notes:
719   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
720   to be used during the `TS` integration.
721 
722 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSSetSolutionFunction()`
723 @*/
724 PetscErrorCode TSMonitorDrawError(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dummy)
725 {
726   TSMonitorDrawCtx ctx    = (TSMonitorDrawCtx)dummy;
727   PetscViewer      viewer = ctx->viewer;
728   Vec              work;
729 
730   PetscFunctionBegin;
731   if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(PETSC_SUCCESS);
732   PetscCall(VecDuplicate(u, &work));
733   PetscCall(TSComputeSolutionFunction(ts, ptime, work));
734   PetscCall(VecAXPY(work, -1.0, u));
735   PetscCall(VecView(work, viewer));
736   PetscCall(VecDestroy(&work));
737   PetscFunctionReturn(PETSC_SUCCESS);
738 }
739 
740 /*@C
741   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
742 
743   Collective
744 
745   Input Parameters:
746 + ts    - the `TS` context
747 . step  - current time-step
748 . ptime - current time
749 . u     - current state
750 - vf    - viewer and its format
751 
752   Level: intermediate
753 
754   Notes:
755   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
756   to be used during the `TS` integration.
757 
758 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`
759 @*/
760 PetscErrorCode TSMonitorSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscViewerAndFormat *vf)
761 {
762   PetscFunctionBegin;
763   if ((vf->view_interval > 0 && !(step % vf->view_interval)) || (vf->view_interval && ts->reason)) {
764     PetscCall(PetscViewerPushFormat(vf->viewer, vf->format));
765     PetscCall(VecView(u, vf->viewer));
766     PetscCall(PetscViewerPopFormat(vf->viewer));
767   }
768   PetscFunctionReturn(PETSC_SUCCESS);
769 }
770 
771 /*@C
772   TSMonitorSolutionVTK - Monitors progress of the `TS` solvers by `VecView()` for the solution at selected timesteps.
773 
774   Collective
775 
776   Input Parameters:
777 + ts    - the `TS` context
778 . step  - current time-step
779 . ptime - current time
780 . u     - current state
781 - ctx   - monitor context obtained with `TSMonitorSolutionVTKCtxCreate()`
782 
783   Level: developer
784 
785   Notes:
786   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.
787   These are named according to the file name template.
788 
789   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
790   to be used during the `TS` integration.
791 
792 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`
793 @*/
794 PetscErrorCode TSMonitorSolutionVTK(TS ts, PetscInt step, PetscReal ptime, Vec u, TSMonitorVTKCtx ctx)
795 {
796   char        filename[PETSC_MAX_PATH_LEN];
797   PetscViewer viewer;
798 
799   PetscFunctionBegin;
800   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
801   if (((ctx->interval > 0) && (!(step % ctx->interval))) || (ctx->interval && ts->reason)) {
802     PetscCall(PetscSNPrintf(filename, sizeof(filename), (const char *)ctx->filenametemplate, step));
803     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts), filename, FILE_MODE_WRITE, &viewer));
804     PetscCall(VecView(u, viewer));
805     PetscCall(PetscViewerDestroy(&viewer));
806   }
807   PetscFunctionReturn(PETSC_SUCCESS);
808 }
809 
810 /*@C
811   TSMonitorSolutionVTKDestroy - Destroy the monitor context created with `TSMonitorSolutionVTKCtxCreate()`
812 
813   Not Collective
814 
815   Input Parameter:
816 . ctx - the monitor context
817 
818   Level: developer
819 
820   Note:
821   This function is normally passed to `TSMonitorSet()` along with `TSMonitorSolutionVTK()`.
822 
823 .seealso: [](ch_ts), `TSMonitorSet()`, `TSMonitorSolutionVTK()`
824 @*/
825 PetscErrorCode TSMonitorSolutionVTKDestroy(TSMonitorVTKCtx *ctx)
826 {
827   PetscFunctionBegin;
828   PetscCall(PetscFree((*ctx)->filenametemplate));
829   PetscCall(PetscFree(*ctx));
830   PetscFunctionReturn(PETSC_SUCCESS);
831 }
832 
833 /*@C
834   TSMonitorSolutionVTKCtxCreate - Create the monitor context to be used in `TSMonitorSolutionVTK()`
835 
836   Not collective
837 
838   Input Parameter:
839 . filenametemplate - the template file name, e.g. foo-%03d.vts
840 
841   Output Parameter:
842 . ctx - the monitor context
843 
844   Level: developer
845 
846   Note:
847   This function is normally used inside `TSSetFromOptions()` to pass the context created to `TSMonitorSet()` along with `TSMonitorSolutionVTK()`.
848 
849 .seealso: [](ch_ts), `TSMonitorSet()`, `TSMonitorSolutionVTK()`, `TSMonitorSolutionVTKDestroy()`
850 @*/
851 PetscErrorCode TSMonitorSolutionVTKCtxCreate(const char *filenametemplate, TSMonitorVTKCtx *ctx)
852 {
853   const char     *ptr = NULL, *ptr2 = NULL;
854   TSMonitorVTKCtx ictx;
855 
856   PetscFunctionBegin;
857   PetscAssertPointer(filenametemplate, 1);
858   PetscAssertPointer(ctx, 2);
859   /* Do some cursory validation of the input. */
860   PetscCall(PetscStrstr(filenametemplate, "%", (char **)&ptr));
861   PetscCheck(ptr, PETSC_COMM_SELF, PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
862   for (ptr++; ptr && *ptr; ptr++) {
863     PetscCall(PetscStrchr("DdiouxX", *ptr, (char **)&ptr2));
864     PetscCheck(ptr2 || (*ptr >= '0' && *ptr <= '9'), PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03" PetscInt_FMT ".vts");
865     if (ptr2) break;
866   }
867   PetscCall(PetscNew(&ictx));
868   PetscCall(PetscStrallocpy(filenametemplate, &ictx->filenametemplate));
869   ictx->interval = 1;
870 
871   *ctx = ictx;
872   PetscFunctionReturn(PETSC_SUCCESS);
873 }
874 
875 /*@C
876   TSMonitorLGSolution - Monitors progress of the `TS` solvers by plotting each component of the solution vector
877   in a time based line graph
878 
879   Collective
880 
881   Input Parameters:
882 + ts    - the `TS` context
883 . step  - current time-step
884 . ptime - current time
885 . u     - current solution
886 - dctx  - the `TSMonitorLGCtx` object that contains all the options for the monitoring, this is created with `TSMonitorLGCtxCreate()`
887 
888   Options Database Key:
889 . -ts_monitor_lg_solution_variables - enable monitor of lg solution variables
890 
891   Level: intermediate
892 
893   Notes:
894   Each process in a parallel run displays its component solutions in a separate window
895 
896   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
897   to be used during the `TS` integration.
898 
899 .seealso: [](ch_ts), `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGCtxCreate()`, `TSMonitorLGCtxSetVariableNames()`, `TSMonitorLGCtxGetVariableNames()`,
900           `TSMonitorLGSetVariableNames()`, `TSMonitorLGGetVariableNames()`, `TSMonitorLGSetDisplayVariables()`, `TSMonitorLGCtxSetDisplayVariables()`,
901           `TSMonitorLGCtxSetTransform()`, `TSMonitorLGSetTransform()`, `TSMonitorLGError()`, `TSMonitorLGSNESIterations()`, `TSMonitorLGKSPIterations()`,
902           `TSMonitorEnvelopeCtxCreate()`, `TSMonitorEnvelopeGetBounds()`, `TSMonitorEnvelopeCtxDestroy()`, `TSMonitorEnvelop()`
903 @*/
904 PetscErrorCode TSMonitorLGSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dctx)
905 {
906   TSMonitorLGCtx     ctx = (TSMonitorLGCtx)dctx;
907   const PetscScalar *yy;
908   Vec                v;
909 
910   PetscFunctionBegin;
911   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
912   if (!step) {
913     PetscDrawAxis axis;
914     PetscInt      dim;
915     PetscCall(PetscDrawLGGetAxis(ctx->lg, &axis));
916     PetscCall(PetscDrawAxisSetLabels(axis, "Solution as function of time", "Time", "Solution"));
917     if (!ctx->names) {
918       PetscBool flg;
919       /* user provides names of variables to plot but no names has been set so assume names are integer values */
920       PetscCall(PetscOptionsHasName(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_monitor_lg_solution_variables", &flg));
921       if (flg) {
922         PetscInt i, n;
923         char   **names;
924         PetscCall(VecGetSize(u, &n));
925         PetscCall(PetscMalloc1(n + 1, &names));
926         for (i = 0; i < n; i++) {
927           PetscCall(PetscMalloc1(5, &names[i]));
928           PetscCall(PetscSNPrintf(names[i], 5, "%" PetscInt_FMT, i));
929         }
930         names[n]   = NULL;
931         ctx->names = names;
932       }
933     }
934     if (ctx->names && !ctx->displaynames) {
935       char    **displaynames;
936       PetscBool flg;
937       PetscCall(VecGetLocalSize(u, &dim));
938       PetscCall(PetscCalloc1(dim + 1, &displaynames));
939       PetscCall(PetscOptionsGetStringArray(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_monitor_lg_solution_variables", displaynames, &dim, &flg));
940       if (flg) PetscCall(TSMonitorLGCtxSetDisplayVariables(ctx, (const char *const *)displaynames));
941       PetscCall(PetscStrArrayDestroy(&displaynames));
942     }
943     if (ctx->displaynames) {
944       PetscCall(PetscDrawLGSetDimension(ctx->lg, ctx->ndisplayvariables));
945       PetscCall(PetscDrawLGSetLegend(ctx->lg, (const char *const *)ctx->displaynames));
946     } else if (ctx->names) {
947       PetscCall(VecGetLocalSize(u, &dim));
948       PetscCall(PetscDrawLGSetDimension(ctx->lg, dim));
949       PetscCall(PetscDrawLGSetLegend(ctx->lg, (const char *const *)ctx->names));
950     } else {
951       PetscCall(VecGetLocalSize(u, &dim));
952       PetscCall(PetscDrawLGSetDimension(ctx->lg, dim));
953     }
954     PetscCall(PetscDrawLGReset(ctx->lg));
955   }
956 
957   if (!ctx->transform) v = u;
958   else PetscCall((*ctx->transform)(ctx->transformctx, u, &v));
959   PetscCall(VecGetArrayRead(v, &yy));
960   if (ctx->displaynames) {
961     PetscInt i;
962     for (i = 0; i < ctx->ndisplayvariables; i++) ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
963     PetscCall(PetscDrawLGAddCommonPoint(ctx->lg, ptime, ctx->displayvalues));
964   } else {
965 #if defined(PETSC_USE_COMPLEX)
966     PetscInt   i, n;
967     PetscReal *yreal;
968     PetscCall(VecGetLocalSize(v, &n));
969     PetscCall(PetscMalloc1(n, &yreal));
970     for (i = 0; i < n; i++) yreal[i] = PetscRealPart(yy[i]);
971     PetscCall(PetscDrawLGAddCommonPoint(ctx->lg, ptime, yreal));
972     PetscCall(PetscFree(yreal));
973 #else
974     PetscCall(PetscDrawLGAddCommonPoint(ctx->lg, ptime, yy));
975 #endif
976   }
977   PetscCall(VecRestoreArrayRead(v, &yy));
978   if (ctx->transform) PetscCall(VecDestroy(&v));
979 
980   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
981     PetscCall(PetscDrawLGDraw(ctx->lg));
982     PetscCall(PetscDrawLGSave(ctx->lg));
983   }
984   PetscFunctionReturn(PETSC_SUCCESS);
985 }
986 
987 /*@C
988   TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
989 
990   Collective
991 
992   Input Parameters:
993 + ts    - the `TS` context
994 - names - the names of the components, final string must be `NULL`
995 
996   Level: intermediate
997 
998   Notes:
999   If the `TS` object does not have a `TSMonitorLGCtx` associated with it then this function is ignored
1000 
1001 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetDisplayVariables()`, `TSMonitorLGCtxSetVariableNames()`
1002 @*/
1003 PetscErrorCode TSMonitorLGSetVariableNames(TS ts, const char *const *names)
1004 {
1005   PetscInt i;
1006 
1007   PetscFunctionBegin;
1008   for (i = 0; i < ts->numbermonitors; i++) {
1009     if (ts->monitor[i] == TSMonitorLGSolution) {
1010       PetscCall(TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i], names));
1011       break;
1012     }
1013   }
1014   PetscFunctionReturn(PETSC_SUCCESS);
1015 }
1016 
1017 /*@C
1018   TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
1019 
1020   Collective
1021 
1022   Input Parameters:
1023 + ctx   - the `TS` context
1024 - names - the names of the components, final string must be `NULL`
1025 
1026   Level: intermediate
1027 
1028 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetDisplayVariables()`, `TSMonitorLGSetVariableNames()`
1029 @*/
1030 PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx, const char *const *names)
1031 {
1032   PetscFunctionBegin;
1033   PetscCall(PetscStrArrayDestroy(&ctx->names));
1034   PetscCall(PetscStrArrayallocpy(names, &ctx->names));
1035   PetscFunctionReturn(PETSC_SUCCESS);
1036 }
1037 
1038 /*@C
1039   TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot
1040 
1041   Collective
1042 
1043   Input Parameter:
1044 . ts - the `TS` context
1045 
1046   Output Parameter:
1047 . names - the names of the components, final string must be `NULL`
1048 
1049   Level: intermediate
1050 
1051   Note:
1052   If the `TS` object does not have a `TSMonitorLGCtx` associated with it then this function is ignored
1053 
1054 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetDisplayVariables()`
1055 @*/
1056 PetscErrorCode TSMonitorLGGetVariableNames(TS ts, const char *const **names)
1057 {
1058   PetscInt i;
1059 
1060   PetscFunctionBegin;
1061   *names = NULL;
1062   for (i = 0; i < ts->numbermonitors; i++) {
1063     if (ts->monitor[i] == TSMonitorLGSolution) {
1064       TSMonitorLGCtx ctx = (TSMonitorLGCtx)ts->monitorcontext[i];
1065       *names             = (const char *const *)ctx->names;
1066       break;
1067     }
1068   }
1069   PetscFunctionReturn(PETSC_SUCCESS);
1070 }
1071 
1072 /*@C
1073   TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor
1074 
1075   Collective
1076 
1077   Input Parameters:
1078 + ctx          - the `TSMonitorLG` context
1079 - displaynames - the names of the components, final string must be `NULL`
1080 
1081   Level: intermediate
1082 
1083 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetVariableNames()`
1084 @*/
1085 PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx, const char *const *displaynames)
1086 {
1087   PetscInt j = 0, k;
1088 
1089   PetscFunctionBegin;
1090   if (!ctx->names) PetscFunctionReturn(PETSC_SUCCESS);
1091   PetscCall(PetscStrArrayDestroy(&ctx->displaynames));
1092   PetscCall(PetscStrArrayallocpy(displaynames, &ctx->displaynames));
1093   while (displaynames[j]) j++;
1094   ctx->ndisplayvariables = j;
1095   PetscCall(PetscMalloc1(ctx->ndisplayvariables, &ctx->displayvariables));
1096   PetscCall(PetscMalloc1(ctx->ndisplayvariables, &ctx->displayvalues));
1097   j = 0;
1098   while (displaynames[j]) {
1099     k = 0;
1100     while (ctx->names[k]) {
1101       PetscBool flg;
1102       PetscCall(PetscStrcmp(displaynames[j], ctx->names[k], &flg));
1103       if (flg) {
1104         ctx->displayvariables[j] = k;
1105         break;
1106       }
1107       k++;
1108     }
1109     j++;
1110   }
1111   PetscFunctionReturn(PETSC_SUCCESS);
1112 }
1113 
1114 /*@C
1115   TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor
1116 
1117   Collective
1118 
1119   Input Parameters:
1120 + ts           - the `TS` context
1121 - displaynames - the names of the components, final string must be `NULL`
1122 
1123   Level: intermediate
1124 
1125   Note:
1126   If the `TS` object does not have a `TSMonitorLGCtx` associated with it then this function is ignored
1127 
1128 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetVariableNames()`
1129 @*/
1130 PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts, const char *const *displaynames)
1131 {
1132   PetscInt i;
1133 
1134   PetscFunctionBegin;
1135   for (i = 0; i < ts->numbermonitors; i++) {
1136     if (ts->monitor[i] == TSMonitorLGSolution) {
1137       PetscCall(TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i], displaynames));
1138       break;
1139     }
1140   }
1141   PetscFunctionReturn(PETSC_SUCCESS);
1142 }
1143 
1144 /*@C
1145   TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed
1146 
1147   Collective
1148 
1149   Input Parameters:
1150 + ts        - the `TS` context
1151 . transform - the transform function
1152 . destroy   - function to destroy the optional context, see `PetscCtxDestroyFn` for its calling sequence
1153 - tctx      - optional context used by transform function
1154 
1155   Level: intermediate
1156 
1157   Note:
1158   If the `TS` object does not have a `TSMonitorLGCtx` associated with it then this function is ignored
1159 
1160 .seealso: [](ch_ts), `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetVariableNames()`, `TSMonitorLGCtxSetTransform()`, `PetscCtxDestroyFn`
1161 @*/
1162 PetscErrorCode TSMonitorLGSetTransform(TS ts, PetscErrorCode (*transform)(void *, Vec, Vec *), PetscCtxDestroyFn *destroy, void *tctx)
1163 {
1164   PetscInt i;
1165 
1166   PetscFunctionBegin;
1167   for (i = 0; i < ts->numbermonitors; i++) {
1168     if (ts->monitor[i] == TSMonitorLGSolution) PetscCall(TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i], transform, destroy, tctx));
1169   }
1170   PetscFunctionReturn(PETSC_SUCCESS);
1171 }
1172 
1173 /*@C
1174   TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed
1175 
1176   Collective
1177 
1178   Input Parameters:
1179 + tctx      - the `TS` context
1180 . transform - the transform function
1181 . destroy   - function to destroy the optional context, see `PetscCtxDestroyFn` for its calling sequence
1182 - ctx       - optional context used by transform function
1183 
1184   Level: intermediate
1185 
1186 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetVariableNames()`, `TSMonitorLGSetTransform()`, `PetscCtxDestroyFn`
1187 @*/
1188 PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx, PetscErrorCode (*transform)(void *, Vec, Vec *), PetscCtxDestroyFn *destroy, void *tctx)
1189 {
1190   PetscFunctionBegin;
1191   ctx->transform        = transform;
1192   ctx->transformdestroy = destroy;
1193   ctx->transformctx     = tctx;
1194   PetscFunctionReturn(PETSC_SUCCESS);
1195 }
1196 
1197 /*@C
1198   TSMonitorLGError - Monitors progress of the `TS` solvers by plotting each component of the error
1199   in a time based line graph
1200 
1201   Collective
1202 
1203   Input Parameters:
1204 + ts    - the `TS` context
1205 . step  - current time-step
1206 . ptime - current time
1207 . u     - current solution
1208 - dummy - `TSMonitorLGCtx` object created with `TSMonitorLGCtxCreate()`
1209 
1210   Options Database Key:
1211 . -ts_monitor_lg_error - create a graphical monitor of error history
1212 
1213   Level: intermediate
1214 
1215   Notes:
1216   Each process in a parallel run displays its component errors in a separate window
1217 
1218   The user must provide the solution using `TSSetSolutionFunction()` to use this monitor.
1219 
1220   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1221   to be used during the TS integration.
1222 
1223 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSSetSolutionFunction()`
1224 @*/
1225 PetscErrorCode TSMonitorLGError(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dummy)
1226 {
1227   TSMonitorLGCtx     ctx = (TSMonitorLGCtx)dummy;
1228   const PetscScalar *yy;
1229   Vec                y;
1230 
1231   PetscFunctionBegin;
1232   if (!step) {
1233     PetscDrawAxis axis;
1234     PetscInt      dim;
1235     PetscCall(PetscDrawLGGetAxis(ctx->lg, &axis));
1236     PetscCall(PetscDrawAxisSetLabels(axis, "Error in solution as function of time", "Time", "Error"));
1237     PetscCall(VecGetLocalSize(u, &dim));
1238     PetscCall(PetscDrawLGSetDimension(ctx->lg, dim));
1239     PetscCall(PetscDrawLGReset(ctx->lg));
1240   }
1241   PetscCall(VecDuplicate(u, &y));
1242   PetscCall(TSComputeSolutionFunction(ts, ptime, y));
1243   PetscCall(VecAXPY(y, -1.0, u));
1244   PetscCall(VecGetArrayRead(y, &yy));
1245 #if defined(PETSC_USE_COMPLEX)
1246   {
1247     PetscReal *yreal;
1248     PetscInt   i, n;
1249     PetscCall(VecGetLocalSize(y, &n));
1250     PetscCall(PetscMalloc1(n, &yreal));
1251     for (i = 0; i < n; i++) yreal[i] = PetscRealPart(yy[i]);
1252     PetscCall(PetscDrawLGAddCommonPoint(ctx->lg, ptime, yreal));
1253     PetscCall(PetscFree(yreal));
1254   }
1255 #else
1256   PetscCall(PetscDrawLGAddCommonPoint(ctx->lg, ptime, yy));
1257 #endif
1258   PetscCall(VecRestoreArrayRead(y, &yy));
1259   PetscCall(VecDestroy(&y));
1260   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
1261     PetscCall(PetscDrawLGDraw(ctx->lg));
1262     PetscCall(PetscDrawLGSave(ctx->lg));
1263   }
1264   PetscFunctionReturn(PETSC_SUCCESS);
1265 }
1266 
1267 /*@C
1268   TSMonitorSPSwarmSolution - Graphically displays phase plots of `DMSWARM` particles on a scatter plot
1269 
1270   Input Parameters:
1271 + ts    - the `TS` context
1272 . step  - current time-step
1273 . ptime - current time
1274 . u     - current solution
1275 - dctx  - the `TSMonitorSPCtx` object that contains all the options for the monitoring, this is created with `TSMonitorSPCtxCreate()`
1276 
1277   Options Database Keys:
1278 + -ts_monitor_sp_swarm <n>                  - Monitor the solution every n steps, or -1 for plotting only the final solution
1279 . -ts_monitor_sp_swarm_retain <n>           - Retain n old points so we can see the history, or -1 for all points
1280 . -ts_monitor_sp_swarm_multi_species <bool> - Color each species differently
1281 - -ts_monitor_sp_swarm_phase <bool>         - Plot in phase space, as opposed to coordinate space
1282 
1283   Level: intermediate
1284 
1285   Notes:
1286   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1287   to be used during the `TS` integration.
1288 
1289 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `DMSWARM`, `TSMonitorSPCtxCreate()`
1290 @*/
1291 PetscErrorCode TSMonitorSPSwarmSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dctx)
1292 {
1293   TSMonitorSPCtx     ctx = (TSMonitorSPCtx)dctx;
1294   PetscDraw          draw;
1295   DM                 dm, cdm;
1296   const PetscScalar *yy;
1297   PetscInt           Np, p, dim = 2, *species;
1298   PetscReal          species_color;
1299 
1300   PetscFunctionBegin;
1301   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
1302   PetscCall(TSGetDM(ts, &dm));
1303   if (!step) {
1304     PetscDrawAxis axis;
1305     PetscReal     dmboxlower[2], dmboxupper[2];
1306 
1307     PetscCall(TSGetDM(ts, &dm));
1308     PetscCall(DMGetDimension(dm, &dim));
1309     PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Monitor only supports two dimensional fields");
1310     PetscCall(DMSwarmGetCellDM(dm, &cdm));
1311     PetscCall(DMGetBoundingBox(cdm, dmboxlower, dmboxupper));
1312     PetscCall(VecGetLocalSize(u, &Np));
1313     Np /= dim * 2;
1314     PetscCall(PetscDrawSPGetAxis(ctx->sp, &axis));
1315     if (ctx->phase) {
1316       PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "X", "V"));
1317       PetscCall(PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], -10, 10));
1318     } else {
1319       PetscCall(PetscDrawAxisSetLabels(axis, "Particles", "X", "Y"));
1320       PetscCall(PetscDrawAxisSetLimits(axis, dmboxlower[0], dmboxupper[0], dmboxlower[1], dmboxupper[1]));
1321     }
1322     PetscCall(PetscDrawAxisSetHoldLimits(axis, PETSC_TRUE));
1323     PetscCall(PetscDrawSPReset(ctx->sp));
1324   }
1325   if (ctx->multispecies) PetscCall(DMSwarmGetField(dm, "species", NULL, NULL, (void **)&species));
1326   PetscCall(VecGetLocalSize(u, &Np));
1327   Np /= dim * 2;
1328   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
1329     PetscCall(PetscDrawSPGetDraw(ctx->sp, &draw));
1330     if ((ctx->retain == 0) || (ctx->retain > 0 && !(step % ctx->retain))) PetscCall(PetscDrawClear(draw));
1331     PetscCall(PetscDrawFlush(draw));
1332     PetscCall(PetscDrawSPReset(ctx->sp));
1333     PetscCall(VecGetArrayRead(u, &yy));
1334     for (p = 0; p < Np; ++p) {
1335       PetscReal x, y;
1336 
1337       if (ctx->phase) {
1338         x = PetscRealPart(yy[p * dim * 2]);
1339         y = PetscRealPart(yy[p * dim * 2 + dim]);
1340       } else {
1341         x = PetscRealPart(yy[p * dim * 2]);
1342         y = PetscRealPart(yy[p * dim * 2 + 1]);
1343       }
1344       if (ctx->multispecies) {
1345         species_color = species[p] + 2;
1346         PetscCall(PetscDrawSPAddPointColorized(ctx->sp, &x, &y, &species_color));
1347       } else {
1348         PetscCall(PetscDrawSPAddPoint(ctx->sp, &x, &y));
1349       }
1350       PetscCall(PetscDrawSPAddPoint(ctx->sp, &x, &y));
1351     }
1352     PetscCall(VecRestoreArrayRead(u, &yy));
1353     PetscCall(PetscDrawSPDraw(ctx->sp, PETSC_FALSE));
1354     PetscCall(PetscDrawSPSave(ctx->sp));
1355     if (ctx->multispecies) PetscCall(DMSwarmRestoreField(dm, "species", NULL, NULL, (void **)&species));
1356   }
1357   PetscFunctionReturn(PETSC_SUCCESS);
1358 }
1359 
1360 /*@C
1361   TSMonitorHGSwarmSolution - Graphically displays histograms of `DMSWARM` particles
1362 
1363   Input Parameters:
1364 + ts    - the `TS` context
1365 . step  - current time-step
1366 . ptime - current time
1367 . u     - current solution
1368 - dctx  - the `TSMonitorSPCtx` object that contains all the options for the monitoring, this is created with `TSMonitorHGCtxCreate()`
1369 
1370   Options Database Keys:
1371 + -ts_monitor_hg_swarm <n>             - Monitor the solution every n steps, or -1 for plotting only the final solution
1372 . -ts_monitor_hg_swarm_species <num>   - Number of species to histogram
1373 . -ts_monitor_hg_swarm_bins <num>      - Number of histogram bins
1374 - -ts_monitor_hg_swarm_velocity <bool> - Plot in velocity space, as opposed to coordinate space
1375 
1376   Level: intermediate
1377 
1378   Note:
1379   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1380   to be used during the `TS` integration.
1381 
1382 .seealso: `TSMonitorSet()`
1383 @*/
1384 PetscErrorCode TSMonitorHGSwarmSolution(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dctx)
1385 {
1386   TSMonitorHGCtx     ctx = (TSMonitorHGCtx)dctx;
1387   PetscDraw          draw;
1388   DM                 sw;
1389   const PetscScalar *yy;
1390   PetscInt          *species;
1391   PetscInt           dim, d = 0, Np, p, Ns, s;
1392 
1393   PetscFunctionBegin;
1394   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
1395   PetscCall(TSGetDM(ts, &sw));
1396   PetscCall(DMGetDimension(sw, &dim));
1397   PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1398   Ns = PetscMin(Ns, ctx->Ns);
1399   PetscCall(VecGetLocalSize(u, &Np));
1400   Np /= dim * 2;
1401   if (!step) {
1402     PetscDrawAxis axis;
1403     char          title[PETSC_MAX_PATH_LEN];
1404 
1405     for (s = 0; s < Ns; ++s) {
1406       PetscCall(PetscDrawHGGetAxis(ctx->hg[s], &axis));
1407       PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Species %" PetscInt_FMT, s));
1408       if (ctx->velocity) PetscCall(PetscDrawAxisSetLabels(axis, title, "V", "N"));
1409       else PetscCall(PetscDrawAxisSetLabels(axis, title, "X", "N"));
1410     }
1411   }
1412   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
1413     PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1414     for (s = 0; s < Ns; ++s) {
1415       PetscCall(PetscDrawHGReset(ctx->hg[s]));
1416       PetscCall(PetscDrawHGGetDraw(ctx->hg[s], &draw));
1417       PetscCall(PetscDrawClear(draw));
1418       PetscCall(PetscDrawFlush(draw));
1419     }
1420     PetscCall(VecGetArrayRead(u, &yy));
1421     for (p = 0; p < Np; ++p) {
1422       const PetscInt s = species[p] < Ns ? species[p] : 0;
1423       PetscReal      v;
1424 
1425       if (ctx->velocity) v = PetscRealPart(yy[p * dim * 2 + d + dim]);
1426       else v = PetscRealPart(yy[p * dim * 2 + d]);
1427       PetscCall(PetscDrawHGAddValue(ctx->hg[s], v));
1428     }
1429     PetscCall(VecRestoreArrayRead(u, &yy));
1430     for (s = 0; s < Ns; ++s) {
1431       PetscCall(PetscDrawHGDraw(ctx->hg[s]));
1432       PetscCall(PetscDrawHGSave(ctx->hg[s]));
1433     }
1434     PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1435   }
1436   PetscFunctionReturn(PETSC_SUCCESS);
1437 }
1438 
1439 /*@C
1440   TSMonitorError - Monitors progress of the `TS` solvers by printing the 2 norm of the error at each timestep
1441 
1442   Collective
1443 
1444   Input Parameters:
1445 + ts    - the `TS` context
1446 . step  - current time-step
1447 . ptime - current time
1448 . u     - current solution
1449 - vf    - unused context
1450 
1451   Options Database Key:
1452 . -ts_monitor_error - create a graphical monitor of error history
1453 
1454   Level: intermediate
1455 
1456   Notes:
1457   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1458   to be used during the `TS` integration.
1459 
1460   The user must provide the solution using `TSSetSolutionFunction()` to use this monitor.
1461 
1462 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSSetSolutionFunction()`
1463 @*/
1464 PetscErrorCode TSMonitorError(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscViewerAndFormat *vf)
1465 {
1466   DM        dm;
1467   PetscDS   ds = NULL;
1468   PetscInt  Nf = -1, f;
1469   PetscBool flg;
1470 
1471   PetscFunctionBegin;
1472   PetscCall(TSGetDM(ts, &dm));
1473   if (dm) PetscCall(DMGetDS(dm, &ds));
1474   if (ds) PetscCall(PetscDSGetNumFields(ds, &Nf));
1475   if (Nf <= 0) {
1476     Vec       y;
1477     PetscReal nrm;
1478 
1479     PetscCall(VecDuplicate(u, &y));
1480     PetscCall(TSComputeSolutionFunction(ts, ptime, y));
1481     PetscCall(VecAXPY(y, -1.0, u));
1482     PetscCall(PetscObjectTypeCompare((PetscObject)vf->viewer, PETSCVIEWERASCII, &flg));
1483     if (flg) {
1484       PetscCall(VecNorm(y, NORM_2, &nrm));
1485       PetscCall(PetscViewerASCIIPrintf(vf->viewer, "2-norm of error %g\n", (double)nrm));
1486     }
1487     PetscCall(PetscObjectTypeCompare((PetscObject)vf->viewer, PETSCVIEWERDRAW, &flg));
1488     if (flg) PetscCall(VecView(y, vf->viewer));
1489     PetscCall(VecDestroy(&y));
1490   } else {
1491     PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
1492     void    **ctxs;
1493     Vec       v;
1494     PetscReal ferrors[1];
1495 
1496     PetscCall(PetscMalloc2(Nf, &exactFuncs, Nf, &ctxs));
1497     for (f = 0; f < Nf; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]));
1498     PetscCall(DMComputeL2FieldDiff(dm, ptime, exactFuncs, ctxs, u, ferrors));
1499     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04" PetscInt_FMT " time = %-8.4g \t L_2 Error: [", step, (double)ptime));
1500     for (f = 0; f < Nf; ++f) {
1501       if (f > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", "));
1502       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%2.3g", (double)ferrors[f]));
1503     }
1504     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "]\n"));
1505 
1506     PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
1507 
1508     PetscCall(PetscOptionsHasName(NULL, NULL, "-exact_vec_view", &flg));
1509     if (flg) {
1510       PetscCall(DMGetGlobalVector(dm, &v));
1511       PetscCall(DMProjectFunction(dm, ptime, exactFuncs, ctxs, INSERT_ALL_VALUES, v));
1512       PetscCall(PetscObjectSetName((PetscObject)v, "Exact Solution"));
1513       PetscCall(VecViewFromOptions(v, NULL, "-exact_vec_view"));
1514       PetscCall(DMRestoreGlobalVector(dm, &v));
1515     }
1516     PetscCall(PetscFree2(exactFuncs, ctxs));
1517   }
1518   PetscFunctionReturn(PETSC_SUCCESS);
1519 }
1520 
1521 PetscErrorCode TSMonitorLGSNESIterations(TS ts, PetscInt n, PetscReal ptime, Vec v, void *monctx)
1522 {
1523   TSMonitorLGCtx ctx = (TSMonitorLGCtx)monctx;
1524   PetscReal      x   = ptime, y;
1525   PetscInt       its;
1526 
1527   PetscFunctionBegin;
1528   if (n < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
1529   if (!n) {
1530     PetscDrawAxis axis;
1531     PetscCall(PetscDrawLGGetAxis(ctx->lg, &axis));
1532     PetscCall(PetscDrawAxisSetLabels(axis, "Nonlinear iterations as function of time", "Time", "SNES Iterations"));
1533     PetscCall(PetscDrawLGReset(ctx->lg));
1534     ctx->snes_its = 0;
1535   }
1536   PetscCall(TSGetSNESIterations(ts, &its));
1537   y = its - ctx->snes_its;
1538   PetscCall(PetscDrawLGAddPoint(ctx->lg, &x, &y));
1539   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
1540     PetscCall(PetscDrawLGDraw(ctx->lg));
1541     PetscCall(PetscDrawLGSave(ctx->lg));
1542   }
1543   ctx->snes_its = its;
1544   PetscFunctionReturn(PETSC_SUCCESS);
1545 }
1546 
1547 PetscErrorCode TSMonitorLGKSPIterations(TS ts, PetscInt n, PetscReal ptime, Vec v, void *monctx)
1548 {
1549   TSMonitorLGCtx ctx = (TSMonitorLGCtx)monctx;
1550   PetscReal      x   = ptime, y;
1551   PetscInt       its;
1552 
1553   PetscFunctionBegin;
1554   if (n < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
1555   if (!n) {
1556     PetscDrawAxis axis;
1557     PetscCall(PetscDrawLGGetAxis(ctx->lg, &axis));
1558     PetscCall(PetscDrawAxisSetLabels(axis, "Linear iterations as function of time", "Time", "KSP Iterations"));
1559     PetscCall(PetscDrawLGReset(ctx->lg));
1560     ctx->ksp_its = 0;
1561   }
1562   PetscCall(TSGetKSPIterations(ts, &its));
1563   y = its - ctx->ksp_its;
1564   PetscCall(PetscDrawLGAddPoint(ctx->lg, &x, &y));
1565   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
1566     PetscCall(PetscDrawLGDraw(ctx->lg));
1567     PetscCall(PetscDrawLGSave(ctx->lg));
1568   }
1569   ctx->ksp_its = its;
1570   PetscFunctionReturn(PETSC_SUCCESS);
1571 }
1572 
1573 /*@C
1574   TSMonitorEnvelopeCtxCreate - Creates a context for use with `TSMonitorEnvelope()`
1575 
1576   Collective
1577 
1578   Input Parameter:
1579 . ts - the `TS` solver object
1580 
1581   Output Parameter:
1582 . ctx - the context
1583 
1584   Level: intermediate
1585 
1586 .seealso: [](ch_ts), `TS`, `TSMonitorLGTimeStep()`, `TSMonitorSet()`, `TSMonitorLGSolution()`, `TSMonitorLGError()`
1587 @*/
1588 PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts, TSMonitorEnvelopeCtx *ctx)
1589 {
1590   PetscFunctionBegin;
1591   PetscCall(PetscNew(ctx));
1592   PetscFunctionReturn(PETSC_SUCCESS);
1593 }
1594 
1595 /*@C
1596   TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution
1597 
1598   Collective
1599 
1600   Input Parameters:
1601 + ts    - the `TS` context
1602 . step  - current time-step
1603 . ptime - current time
1604 . u     - current solution
1605 - dctx  - the envelope context
1606 
1607   Options Database Key:
1608 . -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
1609 
1610   Level: intermediate
1611 
1612   Notes:
1613   After a solve you can use `TSMonitorEnvelopeGetBounds()` to access the envelope
1614 
1615   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1616   to be used during the `TS` integration.
1617 
1618 .seealso: [](ch_ts), `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorEnvelopeGetBounds()`, `TSMonitorEnvelopeCtxCreate()`
1619 @*/
1620 PetscErrorCode TSMonitorEnvelope(TS ts, PetscInt step, PetscReal ptime, Vec u, void *dctx)
1621 {
1622   TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;
1623 
1624   PetscFunctionBegin;
1625   if (!ctx->max) {
1626     PetscCall(VecDuplicate(u, &ctx->max));
1627     PetscCall(VecDuplicate(u, &ctx->min));
1628     PetscCall(VecCopy(u, ctx->max));
1629     PetscCall(VecCopy(u, ctx->min));
1630   } else {
1631     PetscCall(VecPointwiseMax(ctx->max, u, ctx->max));
1632     PetscCall(VecPointwiseMin(ctx->min, u, ctx->min));
1633   }
1634   PetscFunctionReturn(PETSC_SUCCESS);
1635 }
1636 
1637 /*@C
1638   TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution
1639 
1640   Collective
1641 
1642   Input Parameter:
1643 . ts - the `TS` context
1644 
1645   Output Parameters:
1646 + max - the maximum values
1647 - min - the minimum values
1648 
1649   Level: intermediate
1650 
1651   Notes:
1652   If the `TS` does not have a `TSMonitorEnvelopeCtx` associated with it then this function is ignored
1653 
1654 .seealso: [](ch_ts), `TSMonitorEnvelopeCtx`, `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `VecView()`, `TSMonitorLGSetDisplayVariables()`
1655 @*/
1656 PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts, Vec *max, Vec *min)
1657 {
1658   PetscInt i;
1659 
1660   PetscFunctionBegin;
1661   if (max) *max = NULL;
1662   if (min) *min = NULL;
1663   for (i = 0; i < ts->numbermonitors; i++) {
1664     if (ts->monitor[i] == TSMonitorEnvelope) {
1665       TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)ts->monitorcontext[i];
1666       if (max) *max = ctx->max;
1667       if (min) *min = ctx->min;
1668       break;
1669     }
1670   }
1671   PetscFunctionReturn(PETSC_SUCCESS);
1672 }
1673 
1674 /*@C
1675   TSMonitorEnvelopeCtxDestroy - Destroys a context that was created  with `TSMonitorEnvelopeCtxCreate()`.
1676 
1677   Collective
1678 
1679   Input Parameter:
1680 . ctx - the monitor context
1681 
1682   Level: intermediate
1683 
1684 .seealso: [](ch_ts), `TS`, `TSMonitorLGCtxCreate()`, `TSMonitorSet()`, `TSMonitorLGTimeStep()`
1685 @*/
1686 PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
1687 {
1688   PetscFunctionBegin;
1689   PetscCall(VecDestroy(&(*ctx)->min));
1690   PetscCall(VecDestroy(&(*ctx)->max));
1691   PetscCall(PetscFree(*ctx));
1692   PetscFunctionReturn(PETSC_SUCCESS);
1693 }
1694 
1695 /*@C
1696   TSDMSwarmMonitorMoments - Monitors the first three moments of a `DMSWARM` being evolved by the `TS`
1697 
1698   Not Collective
1699 
1700   Input Parameters:
1701 + ts   - the `TS` context
1702 . step - current timestep
1703 . t    - current time
1704 . U    - current solution
1705 - vf   - not used
1706 
1707   Options Database Key:
1708 . -ts_dmswarm_monitor_moments - Monitor moments of particle distribution
1709 
1710   Level: intermediate
1711 
1712   Notes:
1713   This requires a `DMSWARM` be attached to the `TS`.
1714 
1715   This is not called directly by users, rather one calls `TSMonitorSet()`, with this function as an argument, to cause the monitor
1716   to be used during the TS integration.
1717 
1718 .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorDefault()`, `DMSWARM`
1719 @*/
1720 PetscErrorCode TSDMSwarmMonitorMoments(TS ts, PetscInt step, PetscReal t, Vec U, PetscViewerAndFormat *vf)
1721 {
1722   DM                 sw;
1723   const PetscScalar *u;
1724   PetscReal          m = 1.0, totE = 0., totMom[3] = {0., 0., 0.};
1725   PetscInt           dim, d, Np, p;
1726   MPI_Comm           comm;
1727 
1728   PetscFunctionBeginUser;
1729   (void)t;
1730   (void)vf;
1731   PetscCall(TSGetDM(ts, &sw));
1732   if (!sw || step % ts->monitorFrequency != 0) PetscFunctionReturn(PETSC_SUCCESS);
1733   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1734   PetscCall(DMGetDimension(sw, &dim));
1735   PetscCall(VecGetLocalSize(U, &Np));
1736   Np /= dim;
1737   PetscCall(VecGetArrayRead(U, &u));
1738   for (p = 0; p < Np; ++p) {
1739     for (d = 0; d < dim; ++d) {
1740       totE += PetscRealPart(u[p * dim + d] * u[p * dim + d]);
1741       totMom[d] += PetscRealPart(u[p * dim + d]);
1742     }
1743   }
1744   PetscCall(VecRestoreArrayRead(U, &u));
1745   for (d = 0; d < dim; ++d) totMom[d] *= m;
1746   totE *= 0.5 * m;
1747   PetscCall(PetscPrintf(comm, "Step %4" PetscInt_FMT " Total Energy: %10.8lf", step, (double)totE));
1748   for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(comm, "    Total Momentum %c: %10.8lf", (char)('x' + d), (double)totMom[d]));
1749   PetscCall(PetscPrintf(comm, "\n"));
1750   PetscFunctionReturn(PETSC_SUCCESS);
1751 }
1752