Actual source code: tseig.c
petsc-3.9.4 2018-09-11
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: .keywords: TS, monitor, line graph, residual, seealso
46: .seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
48: @*/
49: PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx *ctx)
50: {
51: PetscDraw win;
53: PC pc;
56: PetscNew(ctx);
57: PetscRandomCreate(comm,&(*ctx)->rand);
58: PetscRandomSetFromOptions((*ctx)->rand);
59: PetscDrawCreate(comm,host,label,x,y,m,n,&win);
60: PetscDrawSetFromOptions(win);
61: PetscDrawSPCreate(win,1,&(*ctx)->drawsp);
62: KSPCreate(comm,&(*ctx)->ksp);
63: KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
64: KSPSetType((*ctx)->ksp,KSPGMRES);
65: KSPGMRESSetRestart((*ctx)->ksp,200);
66: KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);
67: KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);
68: KSPSetFromOptions((*ctx)->ksp);
69: KSPGetPC((*ctx)->ksp,&pc);
70: PCSetType(pc,PCNONE);
72: (*ctx)->howoften = howoften;
73: (*ctx)->computeexplicitly = PETSC_FALSE;
75: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);
77: (*ctx)->comm = comm;
78: (*ctx)->xmin = -2.1;
79: (*ctx)->xmax = 1.1;
80: (*ctx)->ymin = -1.1;
81: (*ctx)->ymax = 1.1;
82: return(0);
83: }
85: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg)
86: {
88: PetscReal yr,yi;
91: TSComputeLinearStability(ts,xr,xi,&yr,&yi);
92: if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE;
93: else *flg = PETSC_FALSE;
94: return(0);
95: }
97: PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
98: {
99: TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx;
100: PetscErrorCode ierr;
101: KSP ksp = ctx->ksp;
102: PetscInt n,N,nits,neig,i,its = 200;
103: PetscReal *r,*c,time_step_save;
104: PetscDrawSP drawsp = ctx->drawsp;
105: Mat A,B;
106: Vec xdot;
107: SNES snes;
110: if (step < 0) return(0); /* -1 indicates interpolated solution */
111: if (!step) return(0);
112: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
113: VecDuplicate(v,&xdot);
114: TSGetSNES(ts,&snes);
115: SNESGetJacobian(snes,&A,&B,NULL,NULL);
116: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
117: /*
118: This doesn't work because methods keep and use internal information about the shift so it
119: seems we would need code for each method to trick the correct Jacobian in being computed.
120: */
121: time_step_save = ts->time_step;
122: ts->time_step = PETSC_MAX_REAL;
124: SNESComputeJacobian(snes,v,A,B);
126: ts->time_step = time_step_save;
128: KSPSetOperators(ksp,B,B);
129: VecGetSize(v,&n);
130: if (n < 200) its = n;
131: KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);
132: VecSetRandom(xdot,ctx->rand);
133: KSPSolve(ksp,xdot,xdot);
134: VecDestroy(&xdot);
135: KSPGetIterationNumber(ksp,&nits);
136: N = nits+2;
138: if (nits) {
139: PetscDraw draw;
140: PetscReal pause;
141: PetscDrawAxis axis;
142: PetscReal xmin,xmax,ymin,ymax;
144: PetscDrawSPReset(drawsp);
145: PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);
146: PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c);
147: if (ctx->computeexplicitly) {
148: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
149: neig = n;
150: } else {
151: KSPComputeEigenvalues(ksp,N,r,c,&neig);
152: }
153: /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
154: for (i=0; i<neig; i++) r[i] = -r[i];
155: for (i=0; i<neig; i++) {
156: if (ts->ops->linearstability) {
157: PetscReal fr,fi;
158: TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);
159: if ((fr*fr + fi*fi) > 1.0) {
160: 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));
161: }
162: }
163: PetscDrawSPAddPoint(drawsp,r+i,c+i);
164: }
165: PetscFree2(r,c);
166: PetscDrawSPGetDraw(drawsp,&draw);
167: PetscDrawGetPause(draw,&pause);
168: PetscDrawSetPause(draw,0.0);
169: PetscDrawSPDraw(drawsp,PETSC_TRUE);
170: PetscDrawSetPause(draw,pause);
171: if (ts->ops->linearstability) {
172: PetscDrawSPGetAxis(drawsp,&axis);
173: PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);
174: PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);
175: PetscDrawSPDraw(drawsp,PETSC_FALSE);
176: }
177: PetscDrawSPSave(drawsp);
178: }
179: MatDestroy(&B);
180: }
181: return(0);
182: }
184: /*@C
185: TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate().
187: Collective on TSMonitorSPEigCtx
189: Input Parameter:
190: . ctx - the monitor context
192: Level: intermediate
194: .keywords: TS, monitor, line graph, destroy
196: .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig();
197: @*/
198: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
199: {
200: PetscDraw draw;
204: PetscDrawSPGetDraw((*ctx)->drawsp,&draw);
205: PetscDrawDestroy(&draw);
206: PetscDrawSPDestroy(&(*ctx)->drawsp);
207: KSPDestroy(&(*ctx)->ksp);
208: PetscRandomDestroy(&(*ctx)->rand);
209: PetscFree(*ctx);
210: return(0);
211: }