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