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 on TS 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 Notes: 35 Use TSMonitorSPEigCtxDestroy() to destroy. 36 37 Currently only works if the Jacobian is provided explicitly. 38 39 Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix. 40 41 Level: intermediate 42 43 .seealso: `TSMonitorSPEigTimeStep()`, `TSMonitorSet()`, `TSMonitorLGSolution()`, `TSMonitorLGError()` 44 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 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 PetscReal yr, yi; 82 83 PetscFunctionBegin; 84 PetscCall(TSComputeLinearStability(ts, xr, xi, &yr, &yi)); 85 if ((yr * yr + yi * yi) <= 1.0) *flg = PETSC_TRUE; 86 else *flg = PETSC_FALSE; 87 PetscFunctionReturn(0); 88 } 89 90 PetscErrorCode TSMonitorSPEig(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx) { 91 TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx)monctx; 92 KSP ksp = ctx->ksp; 93 PetscInt n, N, nits, neig, i, its = 200; 94 PetscReal *r, *c, time_step_save; 95 PetscDrawSP drawsp = ctx->drawsp; 96 Mat A, B; 97 Vec xdot; 98 SNES snes; 99 100 PetscFunctionBegin; 101 if (step < 0) PetscFunctionReturn(0); /* -1 indicates interpolated solution */ 102 if (!step) PetscFunctionReturn(0); 103 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 104 PetscCall(VecDuplicate(v, &xdot)); 105 PetscCall(TSGetSNES(ts, &snes)); 106 PetscCall(SNESGetJacobian(snes, &A, &B, NULL, NULL)); 107 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B)); 108 /* 109 This doesn't work because methods keep and use internal information about the shift so it 110 seems we would need code for each method to trick the correct Jacobian in being computed. 111 */ 112 time_step_save = ts->time_step; 113 ts->time_step = PETSC_MAX_REAL; 114 115 PetscCall(SNESComputeJacobian(snes, v, A, B)); 116 117 ts->time_step = time_step_save; 118 119 PetscCall(KSPSetOperators(ksp, B, B)); 120 PetscCall(VecGetSize(v, &n)); 121 if (n < 200) its = n; 122 PetscCall(KSPSetTolerances(ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, its)); 123 PetscCall(VecSetRandom(xdot, ctx->rand)); 124 PetscCall(KSPSolve(ksp, xdot, xdot)); 125 PetscCall(VecDestroy(&xdot)); 126 PetscCall(KSPGetIterationNumber(ksp, &nits)); 127 N = nits + 2; 128 129 if (nits) { 130 PetscDraw draw; 131 PetscReal pause; 132 PetscDrawAxis axis; 133 PetscReal xmin, xmax, ymin, ymax; 134 135 PetscCall(PetscDrawSPReset(drawsp)); 136 PetscCall(PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax)); 137 PetscCall(PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c)); 138 if (ctx->computeexplicitly) { 139 PetscCall(KSPComputeEigenvaluesExplicitly(ksp, n, r, c)); 140 neig = n; 141 } else { 142 PetscCall(KSPComputeEigenvalues(ksp, N, r, c, &neig)); 143 } 144 /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */ 145 for (i = 0; i < neig; i++) r[i] = -r[i]; 146 for (i = 0; i < neig; i++) { 147 if (ts->ops->linearstability) { 148 PetscReal fr, fi; 149 PetscCall(TSComputeLinearStability(ts, r[i], c[i], &fr, &fi)); 150 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))); 151 } 152 PetscCall(PetscDrawSPAddPoint(drawsp, r + i, c + i)); 153 } 154 PetscCall(PetscFree2(r, c)); 155 PetscCall(PetscDrawSPGetDraw(drawsp, &draw)); 156 PetscCall(PetscDrawGetPause(draw, &pause)); 157 PetscCall(PetscDrawSetPause(draw, 0.0)); 158 PetscCall(PetscDrawSPDraw(drawsp, PETSC_TRUE)); 159 PetscCall(PetscDrawSetPause(draw, pause)); 160 if (ts->ops->linearstability) { 161 PetscCall(PetscDrawSPGetAxis(drawsp, &axis)); 162 PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax)); 163 PetscCall(PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode(*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts)); 164 PetscCall(PetscDrawSPDraw(drawsp, PETSC_FALSE)); 165 } 166 PetscCall(PetscDrawSPSave(drawsp)); 167 } 168 PetscCall(MatDestroy(&B)); 169 } 170 PetscFunctionReturn(0); 171 } 172 173 /*@C 174 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 175 176 Collective on TSMonitorSPEigCtx 177 178 Input Parameter: 179 . ctx - the monitor context 180 181 Level: intermediate 182 183 .seealso: `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig();` 184 @*/ 185 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) { 186 PetscDraw draw; 187 188 PetscFunctionBegin; 189 PetscCall(PetscDrawSPGetDraw((*ctx)->drawsp, &draw)); 190 PetscCall(PetscDrawDestroy(&draw)); 191 PetscCall(PetscDrawSPDestroy(&(*ctx)->drawsp)); 192 PetscCall(KSPDestroy(&(*ctx)->ksp)); 193 PetscCall(PetscRandomDestroy(&(*ctx)->rand)); 194 PetscCall(PetscFree(*ctx)); 195 PetscFunctionReturn(0); 196 } 197