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