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