1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdraw.h> 3 4 /* ------------------------------------------------------------------------*/ 5 struct _n_TSMonitorSPEigCtx { 6 PetscDrawSP drawsp; 7 KSP ksp; 8 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 9 PetscBool computeexplicitly; 10 MPI_Comm comm; 11 PetscRandom rand; 12 PetscReal xmin, xmax, ymin, ymax; 13 }; 14 15 /*@C 16 TSMonitorSPEigCtxCreate - Creates a context for use with `TS` to monitor the eigenvalues of the linearized operator 17 18 Collective 19 20 Input Parameters: 21 + comm - the communicator to share the monitor 22 . host - the X display to open, or `NULL` for the local machine 23 . label - the title to put in the title bar 24 . x - the horizontal screen coordinates of the upper left coordinate of the window 25 . y - the vertical coordinates of the upper left coordinate of the window 26 . m - the screen width in pixels 27 . n - the screen height in pixels 28 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 29 30 Output Parameter: 31 . ctx - the context 32 33 Options Database Key: 34 . -ts_monitor_sp_eig - plot egienvalues of linearized right hand side 35 36 Level: intermediate 37 38 Notes: 39 Use `TSMonitorSPEigCtxDestroy()` to destroy the context 40 41 Currently only works if the Jacobian is provided explicitly. 42 43 Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix. 44 45 .seealso: [](ch_ts), `TSMonitorSPEigTimeStep()`, `TSMonitorSet()`, `TSMonitorLGSolution()`, `TSMonitorLGError()` 46 @*/ 47 PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm, const char host[], const char label[], int x, int y, int m, int n, PetscInt howoften, TSMonitorSPEigCtx *ctx) 48 { 49 PetscDraw win; 50 PC pc; 51 52 PetscFunctionBegin; 53 PetscCall(PetscNew(ctx)); 54 PetscCall(PetscRandomCreate(comm, &(*ctx)->rand)); 55 PetscCall(PetscRandomSetFromOptions((*ctx)->rand)); 56 PetscCall(PetscDrawCreate(comm, host, label, x, y, m, n, &win)); 57 PetscCall(PetscDrawSetFromOptions(win)); 58 PetscCall(PetscDrawSPCreate(win, 1, &(*ctx)->drawsp)); 59 PetscCall(KSPCreate(comm, &(*ctx)->ksp)); 60 PetscCall(KSPSetOptionsPrefix((*ctx)->ksp, "ts_monitor_sp_eig_")); /* this is wrong, used use also prefix from the TS */ 61 PetscCall(KSPSetType((*ctx)->ksp, KSPGMRES)); 62 PetscCall(KSPGMRESSetRestart((*ctx)->ksp, 200)); 63 PetscCall(KSPSetTolerances((*ctx)->ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, 200)); 64 PetscCall(KSPSetComputeSingularValues((*ctx)->ksp, PETSC_TRUE)); 65 PetscCall(KSPSetFromOptions((*ctx)->ksp)); 66 PetscCall(KSPGetPC((*ctx)->ksp, &pc)); 67 PetscCall(PCSetType(pc, PCNONE)); 68 69 (*ctx)->howoften = howoften; 70 (*ctx)->computeexplicitly = PETSC_FALSE; 71 72 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_sp_eig_explicitly", &(*ctx)->computeexplicitly, NULL)); 73 74 (*ctx)->comm = comm; 75 (*ctx)->xmin = -2.1; 76 (*ctx)->xmax = 1.1; 77 (*ctx)->ymin = -1.1; 78 (*ctx)->ymax = 1.1; 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr, PetscReal xi, PetscBool *flg) 83 { 84 PetscReal yr, yi; 85 86 PetscFunctionBegin; 87 PetscCall(TSComputeLinearStability(ts, xr, xi, &yr, &yi)); 88 if ((yr * yr + yi * yi) <= 1.0) *flg = PETSC_TRUE; 89 else *flg = PETSC_FALSE; 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 PetscErrorCode TSMonitorSPEig(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx) 94 { 95 TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx)monctx; 96 KSP ksp = ctx->ksp; 97 PetscInt n, N, nits, neig, i, its = 200; 98 PetscReal *r, *c, time_step_save; 99 PetscDrawSP drawsp = ctx->drawsp; 100 Mat A, B; 101 Vec xdot; 102 SNES snes; 103 104 PetscFunctionBegin; 105 if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */ 106 if (!step) PetscFunctionReturn(PETSC_SUCCESS); 107 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 108 PetscCall(VecDuplicate(v, &xdot)); 109 PetscCall(TSGetSNES(ts, &snes)); 110 PetscCall(SNESGetJacobian(snes, &A, &B, NULL, NULL)); 111 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B)); 112 /* 113 This doesn't work because methods keep and use internal information about the shift so it 114 seems we would need code for each method to trick the correct Jacobian in being computed. 115 */ 116 time_step_save = ts->time_step; 117 ts->time_step = PETSC_MAX_REAL; 118 119 PetscCall(SNESComputeJacobian(snes, v, A, B)); 120 121 ts->time_step = time_step_save; 122 123 PetscCall(KSPSetOperators(ksp, B, B)); 124 PetscCall(VecGetSize(v, &n)); 125 if (n < 200) its = n; 126 PetscCall(KSPSetTolerances(ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, its)); 127 PetscCall(VecSetRandom(xdot, ctx->rand)); 128 PetscCall(KSPSolve(ksp, xdot, xdot)); 129 PetscCall(VecDestroy(&xdot)); 130 PetscCall(KSPGetIterationNumber(ksp, &nits)); 131 N = nits + 2; 132 133 if (nits) { 134 PetscDraw draw; 135 PetscReal pause; 136 PetscDrawAxis axis; 137 PetscReal xmin, xmax, ymin, ymax; 138 139 PetscCall(PetscDrawSPReset(drawsp)); 140 PetscCall(PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax)); 141 PetscCall(PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c)); 142 if (ctx->computeexplicitly) { 143 PetscCall(KSPComputeEigenvaluesExplicitly(ksp, n, r, c)); 144 neig = n; 145 } else { 146 PetscCall(KSPComputeEigenvalues(ksp, N, r, c, &neig)); 147 } 148 /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */ 149 for (i = 0; i < neig; i++) r[i] = -r[i]; 150 for (i = 0; i < neig; i++) { 151 if (ts->ops->linearstability) { 152 PetscReal fr, fi; 153 PetscCall(TSComputeLinearStability(ts, r[i], c[i], &fr, &fi)); 154 if ((fr * fr + fi * fi) > 1.0) PetscCall(PetscPrintf(ctx->comm, "Linearized Eigenvalue %g + %g i linear stability function %g norm indicates unstable scheme \n", (double)r[i], (double)c[i], (double)(fr * fr + fi * fi))); 155 } 156 PetscCall(PetscDrawSPAddPoint(drawsp, r + i, c + i)); 157 } 158 PetscCall(PetscFree2(r, c)); 159 PetscCall(PetscDrawSPGetDraw(drawsp, &draw)); 160 PetscCall(PetscDrawGetPause(draw, &pause)); 161 PetscCall(PetscDrawSetPause(draw, 0.0)); 162 PetscCall(PetscDrawSPDraw(drawsp, PETSC_TRUE)); 163 PetscCall(PetscDrawSetPause(draw, pause)); 164 if (ts->ops->linearstability) { 165 PetscCall(PetscDrawSPGetAxis(drawsp, &axis)); 166 PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax)); 167 PetscCall(PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode(*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts)); 168 PetscCall(PetscDrawSPDraw(drawsp, PETSC_FALSE)); 169 } 170 PetscCall(PetscDrawSPSave(drawsp)); 171 } 172 PetscCall(MatDestroy(&B)); 173 } 174 PetscFunctionReturn(PETSC_SUCCESS); 175 } 176 177 /*@C 178 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with `TSMonitorSPEigCtxCreate()`. 179 180 Collective 181 182 Input Parameter: 183 . ctx - the monitor context 184 185 Level: intermediate 186 187 Note: 188 Should be passed to `TSMonitorSet()` along with `TSMonitorSPEig()` an the context created with `TSMonitorSPEigCtxCreate()` 189 190 .seealso: [](ch_ts), `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig()` 191 @*/ 192 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 193 { 194 PetscDraw draw; 195 196 PetscFunctionBegin; 197 PetscCall(PetscDrawSPGetDraw((*ctx)->drawsp, &draw)); 198 PetscCall(PetscDrawDestroy(&draw)); 199 PetscCall(PetscDrawSPDestroy(&(*ctx)->drawsp)); 200 PetscCall(KSPDestroy(&(*ctx)->ksp)); 201 PetscCall(PetscRandomDestroy(&(*ctx)->rand)); 202 PetscCall(PetscFree(*ctx)); 203 PetscFunctionReturn(PETSC_SUCCESS); 204 } 205