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