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