Actual source code: tseig.c
2: #include <petsc/private/tsimpl.h>
3: #include <petscdraw.h>
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: };
16: /*@C
17: TSMonitorSPEigCtxCreate - Creates a context for use with `TS` to monitor the eigenvalues of the linearized operator
19: Collective
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
28: Output Parameter:
29: . ctx - the context
31: Options Database Key:
32: . -ts_monitor_sp_eig - plot egienvalues of linearized right hand side
34: Level: intermediate
36: Notes:
37: Use `TSMonitorSPEigCtxDestroy()` to destroy the context
39: Currently only works if the Jacobian is provided explicitly.
41: Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.
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;
50: PetscNew(ctx);
51: PetscRandomCreate(comm, &(*ctx)->rand);
52: PetscRandomSetFromOptions((*ctx)->rand);
53: PetscDrawCreate(comm, host, label, x, y, m, n, &win);
54: PetscDrawSetFromOptions(win);
55: PetscDrawSPCreate(win, 1, &(*ctx)->drawsp);
56: KSPCreate(comm, &(*ctx)->ksp);
57: KSPSetOptionsPrefix((*ctx)->ksp, "ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
58: KSPSetType((*ctx)->ksp, KSPGMRES);
59: KSPGMRESSetRestart((*ctx)->ksp, 200);
60: KSPSetTolerances((*ctx)->ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, 200);
61: KSPSetComputeSingularValues((*ctx)->ksp, PETSC_TRUE);
62: KSPSetFromOptions((*ctx)->ksp);
63: KSPGetPC((*ctx)->ksp, &pc);
64: PCSetType(pc, PCNONE);
66: (*ctx)->howoften = howoften;
67: (*ctx)->computeexplicitly = PETSC_FALSE;
69: PetscOptionsGetBool(NULL, NULL, "-ts_monitor_sp_eig_explicitly", &(*ctx)->computeexplicitly, NULL);
71: (*ctx)->comm = comm;
72: (*ctx)->xmin = -2.1;
73: (*ctx)->xmax = 1.1;
74: (*ctx)->ymin = -1.1;
75: (*ctx)->ymax = 1.1;
76: return 0;
77: }
79: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr, PetscReal xi, PetscBool *flg)
80: {
81: PetscReal yr, yi;
83: TSComputeLinearStability(ts, xr, xi, &yr, &yi);
84: if ((yr * yr + yi * yi) <= 1.0) *flg = PETSC_TRUE;
85: else *flg = PETSC_FALSE;
86: return 0;
87: }
89: PetscErrorCode TSMonitorSPEig(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx)
90: {
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;
100: if (step < 0) return 0; /* -1 indicates interpolated solution */
101: if (!step) return 0;
102: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
103: VecDuplicate(v, &xdot);
104: TSGetSNES(ts, &snes);
105: SNESGetJacobian(snes, &A, &B, NULL, NULL);
106: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B);
107: /*
108: This doesn't work because methods keep and use internal information about the shift so it
109: seems we would need code for each method to trick the correct Jacobian in being computed.
110: */
111: time_step_save = ts->time_step;
112: ts->time_step = PETSC_MAX_REAL;
114: SNESComputeJacobian(snes, v, A, B);
116: ts->time_step = time_step_save;
118: KSPSetOperators(ksp, B, B);
119: VecGetSize(v, &n);
120: if (n < 200) its = n;
121: KSPSetTolerances(ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, its);
122: VecSetRandom(xdot, ctx->rand);
123: KSPSolve(ksp, xdot, xdot);
124: VecDestroy(&xdot);
125: KSPGetIterationNumber(ksp, &nits);
126: N = nits + 2;
128: if (nits) {
129: PetscDraw draw;
130: PetscReal pause;
131: PetscDrawAxis axis;
132: PetscReal xmin, xmax, ymin, ymax;
134: PetscDrawSPReset(drawsp);
135: PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax);
136: PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c);
137: if (ctx->computeexplicitly) {
138: KSPComputeEigenvaluesExplicitly(ksp, n, r, c);
139: neig = n;
140: } else {
141: KSPComputeEigenvalues(ksp, N, r, c, &neig);
142: }
143: /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
144: for (i = 0; i < neig; i++) r[i] = -r[i];
145: for (i = 0; i < neig; i++) {
146: if (ts->ops->linearstability) {
147: PetscReal fr, fi;
148: TSComputeLinearStability(ts, r[i], c[i], &fr, &fi);
149: if ((fr * fr + fi * fi) > 1.0) 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));
150: }
151: PetscDrawSPAddPoint(drawsp, r + i, c + i);
152: }
153: PetscFree2(r, c);
154: PetscDrawSPGetDraw(drawsp, &draw);
155: PetscDrawGetPause(draw, &pause);
156: PetscDrawSetPause(draw, 0.0);
157: PetscDrawSPDraw(drawsp, PETSC_TRUE);
158: PetscDrawSetPause(draw, pause);
159: if (ts->ops->linearstability) {
160: PetscDrawSPGetAxis(drawsp, &axis);
161: PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax);
162: PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode(*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts);
163: PetscDrawSPDraw(drawsp, PETSC_FALSE);
164: }
165: PetscDrawSPSave(drawsp);
166: }
167: MatDestroy(&B);
168: }
169: return 0;
170: }
172: /*@C
173: TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with `TSMonitorSPEigCtxCreate()`.
175: Collective
177: Input Parameter:
178: . ctx - the monitor context
180: Level: intermediate
182: Note:
183: Should be passed to `TSMonitorSet()` along with `TSMonitorSPEig()` an the context created with `TSMonitorSPEigCtxCreate()`
185: .seealso: [](chapter_ts), `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig();`
186: @*/
187: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
188: {
189: PetscDraw draw;
191: PetscDrawSPGetDraw((*ctx)->drawsp, &draw);
192: PetscDrawDestroy(&draw);
193: PetscDrawSPDestroy(&(*ctx)->drawsp);
194: KSPDestroy(&(*ctx)->ksp);
195: PetscRandomDestroy(&(*ctx)->rand);
196: PetscFree(*ctx);
197: return 0;
198: }