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 @*/
TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx * ctx)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
TSLinearStabilityIndicator(TS ts,PetscReal xr,PetscReal xi,PetscBool * flg)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
TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void * monctx)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 @*/
TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx * ctx)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