Actual source code: tseig.c
petsc-3.14.6 2021-03-30
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: };
17: /*@C
18: TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator
20: Collective on TS
22: Input Parameters:
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, y - the screen coordinates of the upper left coordinate of the window
26: . m, n - the screen width and 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: Notes:
36: Use TSMonitorSPEigCtxDestroy() to destroy.
38: Currently only works if the Jacobian is provided explicitly.
40: Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.
42: Level: intermediate
44: .seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
46: @*/
47: PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx *ctx)
48: {
49: PetscDraw win;
51: PC pc;
54: PetscNew(ctx);
55: PetscRandomCreate(comm,&(*ctx)->rand);
56: PetscRandomSetFromOptions((*ctx)->rand);
57: PetscDrawCreate(comm,host,label,x,y,m,n,&win);
58: PetscDrawSetFromOptions(win);
59: PetscDrawSPCreate(win,1,&(*ctx)->drawsp);
60: KSPCreate(comm,&(*ctx)->ksp);
61: KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
62: KSPSetType((*ctx)->ksp,KSPGMRES);
63: KSPGMRESSetRestart((*ctx)->ksp,200);
64: KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);
65: KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);
66: KSPSetFromOptions((*ctx)->ksp);
67: KSPGetPC((*ctx)->ksp,&pc);
68: PCSetType(pc,PCNONE);
70: (*ctx)->howoften = howoften;
71: (*ctx)->computeexplicitly = PETSC_FALSE;
73: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);
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: return(0);
81: }
83: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg)
84: {
86: PetscReal yr,yi;
89: TSComputeLinearStability(ts,xr,xi,&yr,&yi);
90: if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE;
91: else *flg = PETSC_FALSE;
92: return(0);
93: }
95: PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
96: {
97: TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx;
98: PetscErrorCode ierr;
99: KSP ksp = ctx->ksp;
100: PetscInt n,N,nits,neig,i,its = 200;
101: PetscReal *r,*c,time_step_save;
102: PetscDrawSP drawsp = ctx->drawsp;
103: Mat A,B;
104: Vec xdot;
105: SNES snes;
108: if (step < 0) return(0); /* -1 indicates interpolated solution */
109: if (!step) return(0);
110: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
111: VecDuplicate(v,&xdot);
112: TSGetSNES(ts,&snes);
113: SNESGetJacobian(snes,&A,&B,NULL,NULL);
114: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
115: /*
116: This doesn't work because methods keep and use internal information about the shift so it
117: seems we would need code for each method to trick the correct Jacobian in being computed.
118: */
119: time_step_save = ts->time_step;
120: ts->time_step = PETSC_MAX_REAL;
122: SNESComputeJacobian(snes,v,A,B);
124: ts->time_step = time_step_save;
126: KSPSetOperators(ksp,B,B);
127: VecGetSize(v,&n);
128: if (n < 200) its = n;
129: KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);
130: VecSetRandom(xdot,ctx->rand);
131: KSPSolve(ksp,xdot,xdot);
132: VecDestroy(&xdot);
133: KSPGetIterationNumber(ksp,&nits);
134: N = nits+2;
136: if (nits) {
137: PetscDraw draw;
138: PetscReal pause;
139: PetscDrawAxis axis;
140: PetscReal xmin,xmax,ymin,ymax;
142: PetscDrawSPReset(drawsp);
143: PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);
144: PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c);
145: if (ctx->computeexplicitly) {
146: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
147: neig = n;
148: } else {
149: KSPComputeEigenvalues(ksp,N,r,c,&neig);
150: }
151: /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
152: for (i=0; i<neig; i++) r[i] = -r[i];
153: for (i=0; i<neig; i++) {
154: if (ts->ops->linearstability) {
155: PetscReal fr,fi;
156: TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);
157: if ((fr*fr + fi*fi) > 1.0) {
158: 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));
159: }
160: }
161: PetscDrawSPAddPoint(drawsp,r+i,c+i);
162: }
163: PetscFree2(r,c);
164: PetscDrawSPGetDraw(drawsp,&draw);
165: PetscDrawGetPause(draw,&pause);
166: PetscDrawSetPause(draw,0.0);
167: PetscDrawSPDraw(drawsp,PETSC_TRUE);
168: PetscDrawSetPause(draw,pause);
169: if (ts->ops->linearstability) {
170: PetscDrawSPGetAxis(drawsp,&axis);
171: PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);
172: PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);
173: PetscDrawSPDraw(drawsp,PETSC_FALSE);
174: }
175: PetscDrawSPSave(drawsp);
176: }
177: MatDestroy(&B);
178: }
179: return(0);
180: }
182: /*@C
183: TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate().
185: Collective on TSMonitorSPEigCtx
187: Input Parameter:
188: . ctx - the monitor context
190: Level: intermediate
192: .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig();
193: @*/
194: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
195: {
196: PetscDraw draw;
200: PetscDrawSPGetDraw((*ctx)->drawsp,&draw);
201: PetscDrawDestroy(&draw);
202: PetscDrawSPDestroy(&(*ctx)->drawsp);
203: KSPDestroy(&(*ctx)->ksp);
204: PetscRandomDestroy(&(*ctx)->rand);
205: PetscFree(*ctx);
206: return(0);
207: }