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