Actual source code: tseig.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdraw.h>
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: };
14: /*@C
15: TSMonitorSPEigCtxCreate - Creates a context for use with `TS` to monitor the eigenvalues of the linearized operator
17: Collective
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
29: Output Parameter:
30: . ctx - the context
32: Options Database Key:
33: . -ts_monitor_sp_eig - plot egienvalues of linearized right-hand side
35: Level: intermediate
37: Notes:
38: Use `TSMonitorSPEigCtxDestroy()` to destroy the context
40: Currently only works if the Jacobian is provided explicitly.
42: Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.
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;
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));
68: (*ctx)->howoften = howoften;
69: (*ctx)->computeexplicitly = PETSC_FALSE;
71: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_sp_eig_explicitly", &(*ctx)->computeexplicitly, NULL));
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: }
81: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr, PetscReal xi, PetscBool *flg)
82: {
83: PetscReal yr, yi;
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: }
92: /*@C
93: TSMonitorSPEig - Monitors the eigenvalues of the linearized right-hand-side operator on a scatter plot at each time step
95: Collective
97: Input Parameters:
98: + ts - the `TS` context
99: . step - the current time-step number (a negative value indicates an interpolated solution and is ignored)
100: . ptime - the current time
101: . v - the current solution
102: - monctx - the `TSMonitorSPEigCtx` context, created with `TSMonitorSPEigCtxCreate()`
104: Options Database Key:
105: . -ts_monitor_sp_eig - plot eigenvalues of linearized right-hand side
107: Level: intermediate
109: Notes:
110: This is not called directly by users; pass this function to `TSMonitorSet()` along with the context created by `TSMonitorSPEigCtxCreate()` and `TSMonitorSPEigCtxDestroy()`.
112: Currently only works when the Jacobian is provided explicitly and the ODE has no mass matrix.
114: .seealso: [](ch_ts), `TS`, `TSMonitorSet()`, `TSMonitorSPEigCtxCreate()`, `TSMonitorSPEigCtxDestroy()`
115: @*/
116: PetscErrorCode TSMonitorSPEig(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx)
117: {
118: TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx)monctx;
119: KSP ksp = ctx->ksp;
120: PetscInt n, N, nits, neig, i, its = 200;
121: PetscReal *r, *c, time_step_save;
122: PetscDrawSP drawsp = ctx->drawsp;
123: Mat A, B;
124: Vec xdot;
125: SNES snes;
127: PetscFunctionBegin;
128: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* -1 indicates interpolated solution */
129: if (!step) PetscFunctionReturn(PETSC_SUCCESS);
130: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
131: PetscCall(VecDuplicate(v, &xdot));
132: PetscCall(TSGetSNES(ts, &snes));
133: PetscCall(SNESGetJacobian(snes, &A, &B, NULL, NULL));
134: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
135: /*
136: This doesn't work because methods keep and use internal information about the shift so it
137: seems we would need code for each method to trick the correct Jacobian in being computed.
138: */
139: time_step_save = ts->time_step;
140: ts->time_step = PETSC_MAX_REAL;
142: PetscCall(SNESComputeJacobian(snes, v, A, B));
144: ts->time_step = time_step_save;
146: PetscCall(KSPSetOperators(ksp, B, B));
147: PetscCall(VecGetSize(v, &n));
148: if (n < 200) its = n;
149: PetscCall(KSPSetTolerances(ksp, 1.e-10, PETSC_CURRENT, PETSC_CURRENT, its));
150: PetscCall(VecSetRandom(xdot, ctx->rand));
151: PetscCall(KSPSolve(ksp, xdot, xdot));
152: PetscCall(VecDestroy(&xdot));
153: PetscCall(KSPGetIterationNumber(ksp, &nits));
154: N = nits + 2;
156: if (nits) {
157: PetscDraw draw;
158: PetscReal pause;
159: PetscDrawAxis axis;
160: PetscReal xmin, xmax, ymin, ymax;
162: PetscCall(PetscDrawSPReset(drawsp));
163: PetscCall(PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax));
164: PetscCall(PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c));
165: if (ctx->computeexplicitly) {
166: PetscCall(KSPComputeEigenvaluesExplicitly(ksp, n, r, c));
167: neig = n;
168: } else {
169: PetscCall(KSPComputeEigenvalues(ksp, N, r, c, &neig));
170: }
171: /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
172: for (i = 0; i < neig; i++) r[i] = -r[i];
173: for (i = 0; i < neig; i++) {
174: if (ts->ops->linearstability) {
175: PetscReal fr, fi;
176: PetscCall(TSComputeLinearStability(ts, r[i], c[i], &fr, &fi));
177: 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)));
178: }
179: PetscCall(PetscDrawSPAddPoint(drawsp, r + i, c + i));
180: }
181: PetscCall(PetscFree2(r, c));
182: PetscCall(PetscDrawSPGetDraw(drawsp, &draw));
183: PetscCall(PetscDrawGetPause(draw, &pause));
184: PetscCall(PetscDrawSetPause(draw, 0.0));
185: PetscCall(PetscDrawSPDraw(drawsp, PETSC_TRUE));
186: PetscCall(PetscDrawSetPause(draw, pause));
187: if (ts->ops->linearstability) {
188: PetscCall(PetscDrawSPGetAxis(drawsp, &axis));
189: PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax));
190: PetscCall(PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode (*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts));
191: PetscCall(PetscDrawSPDraw(drawsp, PETSC_FALSE));
192: }
193: PetscCall(PetscDrawSPSave(drawsp));
194: }
195: PetscCall(MatDestroy(&B));
196: }
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*@C
201: TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with `TSMonitorSPEigCtxCreate()`.
203: Collective
205: Input Parameter:
206: . ctx - the monitor context
208: Level: intermediate
210: Note:
211: Should be passed to `TSMonitorSet()` along with `TSMonitorSPEig()` an the context created with `TSMonitorSPEigCtxCreate()`
213: .seealso: [](ch_ts), `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig()`
214: @*/
215: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
216: {
217: PetscDraw draw;
219: PetscFunctionBegin;
220: PetscCall(PetscDrawSPGetDraw((*ctx)->drawsp, &draw));
221: PetscCall(PetscDrawDestroy(&draw));
222: PetscCall(PetscDrawSPDestroy(&(*ctx)->drawsp));
223: PetscCall(KSPDestroy(&(*ctx)->ksp));
224: PetscCall(PetscRandomDestroy(&(*ctx)->rand));
225: PetscCall(PetscFree(*ctx));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }