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