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