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