Actual source code: tseig.c
petsc-3.4.5 2014-06-29
1: #define PETSC_DESIRE_COMPLEX
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
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: };
19: /*@C
20: TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator
22: Collective on TS
24: Input Parameters:
25: + host - the X display to open, or null for the local machine
26: . label - the title to put in the title bar
27: . x, y - the screen coordinates of the upper left coordinate of the window
28: . m, n - the screen width and height in pixels
29: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
31: Output Parameter:
32: . ctx - the context
34: Options Database Key:
35: . -ts_monitor_sp_eig - plot egienvalues of linearized right hand side
37: Notes:
38: Use TSMonitorSPEigCtxDestroy() to destroy.
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: Level: intermediate
46: .keywords: TS, monitor, line graph, residual, seealso
48: .seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
50: @*/
51: PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx *ctx)
52: {
53: PetscDraw win;
55: PC pc;
58: PetscNew(struct _n_TSMonitorSPEigCtx,ctx);
59: PetscRandomCreate(comm,&(*ctx)->rand);
60: PetscRandomSetFromOptions((*ctx)->rand);
61: PetscDrawCreate(comm,host,label,x,y,m,n,&win);
62: PetscDrawSetFromOptions(win);
63: PetscDrawSPCreate(win,1,&(*ctx)->drawsp);
64: KSPCreate(comm,&(*ctx)->ksp);
65: KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
66: KSPSetType((*ctx)->ksp,KSPGMRES);
67: KSPGMRESSetRestart((*ctx)->ksp,200);
68: KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);
69: KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);
70: KSPSetFromOptions((*ctx)->ksp);
71: KSPGetPC((*ctx)->ksp,&pc);
72: PCSetType(pc,PCNONE);
74: (*ctx)->howoften = howoften;
75: (*ctx)->computeexplicitly = PETSC_FALSE;
77: PetscOptionsGetBool(NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);
79: (*ctx)->comm = comm;
80: (*ctx)->xmin = -2.1;
81: (*ctx)->xmax = 1.1;
82: (*ctx)->ymin = -1.1;
83: (*ctx)->ymax = 1.1;
84: return(0);
85: }
89: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg)
90: {
92: PetscReal yr,yi;
95: TSComputeLinearStability(ts,xr,xi,&yr,&yi);
96: if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE;
97: else *flg = PETSC_FALSE;
98: return(0);
99: }
103: PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
104: {
105: TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx;
106: PetscErrorCode ierr;
107: KSP ksp = ctx->ksp;
108: PetscInt n,N,nits,neig,i,its = 200;
109: PetscReal *r,*c,time_step_save;
110: PetscDrawSP drawsp = ctx->drawsp;
111: MatStructure structure;
112: Mat A,B;
113: Vec xdot;
114: SNES snes;
117: if (!step) return(0);
118: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
119: VecDuplicate(v,&xdot);
120: TSGetSNES(ts,&snes);
121: SNESGetJacobian(snes,&A,&B,NULL,NULL);
122: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
123: /*
124: This doesn't work because methods keep and use internal information about the shift so it
125: seems we would need code for each method to trick the correct Jacobian in being computed.
126: */
127: time_step_save = ts->time_step;
128: ts->time_step = PETSC_MAX_REAL;
130: SNESComputeJacobian(snes,v,&A,&B,&structure);
132: ts->time_step = time_step_save;
134: KSPSetOperators(ksp,B,B,structure);
135: VecGetSize(v,&n);
136: if (n < 200) its = n;
137: KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);
138: VecSetRandom(xdot,ctx->rand);
139: KSPSolve(ksp,xdot,xdot);
140: VecDestroy(&xdot);
141: KSPGetIterationNumber(ksp,&nits);
142: N = nits+2;
144: if (nits) {
145: PetscDraw draw;
146: PetscReal pause;
147: PetscDrawAxis axis;
148: PetscReal xmin,xmax,ymin,ymax;
150: PetscDrawSPReset(drawsp);
151: PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);
152: PetscMalloc2(PetscMax(n,N),PetscReal,&r,PetscMax(n,N),PetscReal,&c);
153: if (ctx->computeexplicitly) {
154: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
155: neig = n;
156: } else {
157: KSPComputeEigenvalues(ksp,N,r,c,&neig);
158: }
159: /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
160: for (i=0; i<neig; i++) r[i] = -r[i];
161: for (i=0; i<neig; i++) {
162: if (ts->ops->linearstability) {
163: PetscReal fr,fi;
164: TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);
165: if ((fr*fr + fi*fi) > 1.0) {
166: 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));
167: }
168: }
169: PetscDrawSPAddPoint(drawsp,r+i,c+i);
170: }
171: PetscFree2(r,c);
172: PetscDrawSPGetDraw(drawsp,&draw);
173: PetscDrawGetPause(draw,&pause);
174: PetscDrawSetPause(draw,0.0);
175: PetscDrawSPDraw(drawsp,PETSC_TRUE);
176: PetscDrawSetPause(draw,pause);
178: if (ts->ops->linearstability) {
179: PetscDrawSPGetAxis(drawsp,&axis);
180: PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);
181: PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);
182: PetscDrawSPDraw(drawsp,PETSC_FALSE);
183: }
184: }
185: MatDestroy(&B);
186: }
187: return(0);
188: }
192: /*@C
193: TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate().
195: Collective on TSMonitorSPEigCtx
197: Input Parameter:
198: . ctx - the monitor context
200: Level: intermediate
202: .keywords: TS, monitor, line graph, destroy
204: .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig();
205: @*/
206: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
207: {
208: PetscDraw draw;
212: PetscDrawSPGetDraw((*ctx)->drawsp,&draw);
213: PetscDrawDestroy(&draw);
214: PetscDrawSPDestroy(&(*ctx)->drawsp);
215: KSPDestroy(&(*ctx)->ksp);
216: PetscRandomDestroy(&(*ctx)->rand);
217: PetscFree(*ctx);
218: return(0);
219: }