Actual source code: tseig.c


  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: };

 16: /*@C
 17:    TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator

 19:    Collective on TS

 21:    Input Parameters:
 22: +  host - the X display to open, or null for the local machine
 23: .  label - the title to put in the title bar
 24: .  x, y - the screen coordinates of the upper left coordinate of the window
 25: .  m, n - the screen width and height in pixels
 26: -  howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

 28:    Output Parameter:
 29: .  ctx - the context

 31:    Options Database Key:
 32: .  -ts_monitor_sp_eig - plot egienvalues of linearized right hand side

 34:    Notes:
 35:    Use TSMonitorSPEigCtxDestroy() to destroy.

 37:    Currently only works if the Jacobian is provided explicitly.

 39:    Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.

 41:    Level: intermediate

 43: .seealso: 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;
 50:   PC             pc;

 53:   PetscNew(ctx);
 54:   PetscRandomCreate(comm,&(*ctx)->rand);
 55:   PetscRandomSetFromOptions((*ctx)->rand);
 56:   PetscDrawCreate(comm,host,label,x,y,m,n,&win);
 57:   PetscDrawSetFromOptions(win);
 58:   PetscDrawSPCreate(win,1,&(*ctx)->drawsp);
 59:   KSPCreate(comm,&(*ctx)->ksp);
 60:   KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
 61:   KSPSetType((*ctx)->ksp,KSPGMRES);
 62:   KSPGMRESSetRestart((*ctx)->ksp,200);
 63:   KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);
 64:   KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);
 65:   KSPSetFromOptions((*ctx)->ksp);
 66:   KSPGetPC((*ctx)->ksp,&pc);
 67:   PCSetType(pc,PCNONE);

 69:   (*ctx)->howoften          = howoften;
 70:   (*ctx)->computeexplicitly = PETSC_FALSE;

 72:   PetscOptionsGetBool(NULL,NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);

 74:   (*ctx)->comm = comm;
 75:   (*ctx)->xmin = -2.1;
 76:   (*ctx)->xmax = 1.1;
 77:   (*ctx)->ymin = -1.1;
 78:   (*ctx)->ymax = 1.1;
 79:   return(0);
 80: }

 82: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg)
 83: {
 85:   PetscReal      yr,yi;

 88:   TSComputeLinearStability(ts,xr,xi,&yr,&yi);
 89:   if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE;
 90:   else *flg = PETSC_FALSE;
 91:   return(0);
 92: }

 94: PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
 95: {
 96:   TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx;
 97:   PetscErrorCode    ierr;
 98:   KSP               ksp = ctx->ksp;
 99:   PetscInt          n,N,nits,neig,i,its = 200;
100:   PetscReal         *r,*c,time_step_save;
101:   PetscDrawSP       drawsp = ctx->drawsp;
102:   Mat               A,B;
103:   Vec               xdot;
104:   SNES              snes;

107:   if (step < 0) return(0); /* -1 indicates interpolated solution */
108:   if (!step) return(0);
109:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
110:     VecDuplicate(v,&xdot);
111:     TSGetSNES(ts,&snes);
112:     SNESGetJacobian(snes,&A,&B,NULL,NULL);
113:     MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
114:     /*
115:        This doesn't work because methods keep and use internal information about the shift so it
116:        seems we would need code for each method to trick the correct Jacobian in being computed.
117:      */
118:     time_step_save = ts->time_step;
119:     ts->time_step  = PETSC_MAX_REAL;

121:     SNESComputeJacobian(snes,v,A,B);

123:     ts->time_step  = time_step_save;

125:     KSPSetOperators(ksp,B,B);
126:     VecGetSize(v,&n);
127:     if (n < 200) its = n;
128:     KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);
129:     VecSetRandom(xdot,ctx->rand);
130:     KSPSolve(ksp,xdot,xdot);
131:     VecDestroy(&xdot);
132:     KSPGetIterationNumber(ksp,&nits);
133:     N    = nits+2;

135:     if (nits) {
136:       PetscDraw     draw;
137:       PetscReal     pause;
138:       PetscDrawAxis axis;
139:       PetscReal     xmin,xmax,ymin,ymax;

141:       PetscDrawSPReset(drawsp);
142:       PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);
143:       PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c);
144:       if (ctx->computeexplicitly) {
145:         KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
146:         neig = n;
147:       } else {
148:         KSPComputeEigenvalues(ksp,N,r,c,&neig);
149:       }
150:       /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
151:       for (i=0; i<neig; i++) r[i] = -r[i];
152:       for (i=0; i<neig; i++) {
153:         if (ts->ops->linearstability) {
154:           PetscReal fr,fi;
155:           TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);
156:           if ((fr*fr + fi*fi) > 1.0) {
157:             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));
158:           }
159:         }
160:         PetscDrawSPAddPoint(drawsp,r+i,c+i);
161:       }
162:       PetscFree2(r,c);
163:       PetscDrawSPGetDraw(drawsp,&draw);
164:       PetscDrawGetPause(draw,&pause);
165:       PetscDrawSetPause(draw,0.0);
166:       PetscDrawSPDraw(drawsp,PETSC_TRUE);
167:       PetscDrawSetPause(draw,pause);
168:       if (ts->ops->linearstability) {
169:         PetscDrawSPGetAxis(drawsp,&axis);
170:         PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);
171:         PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);
172:         PetscDrawSPDraw(drawsp,PETSC_FALSE);
173:       }
174:       PetscDrawSPSave(drawsp);
175:     }
176:     MatDestroy(&B);
177:   }
178:   return(0);
179: }

181: /*@C
182:    TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate().

184:    Collective on TSMonitorSPEigCtx

186:    Input Parameter:
187: .  ctx - the monitor context

189:    Level: intermediate

191: .seealso: TSMonitorSPEigCtxCreate(),  TSMonitorSet(), TSMonitorSPEig();
192: @*/
193: PetscErrorCode  TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
194: {
195:   PetscDraw      draw;

199:   PetscDrawSPGetDraw((*ctx)->drawsp,&draw);
200:   PetscDrawDestroy(&draw);
201:   PetscDrawSPDestroy(&(*ctx)->drawsp);
202:   KSPDestroy(&(*ctx)->ksp);
203:   PetscRandomDestroy(&(*ctx)->rand);
204:   PetscFree(*ctx);
205:   return(0);
206: }