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