1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscdraw.h> 4 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 }; 15 16 /*@C 17 TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator 18 19 Collective on TS 20 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 27 28 Output Parameter: 29 . ctx - the context 30 31 Options Database Key: 32 . -ts_monitor_sp_eig - plot egienvalues of linearized right hand side 33 34 Notes: 35 Use TSMonitorSPEigCtxDestroy() to destroy. 36 37 Currently only works if the Jacobian is provided explicitly. 38 39 Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix. 40 41 Level: intermediate 42 43 .seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 44 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; 49 PC pc; 50 51 PetscFunctionBegin; 52 CHKERRQ(PetscNew(ctx)); 53 CHKERRQ(PetscRandomCreate(comm,&(*ctx)->rand)); 54 CHKERRQ(PetscRandomSetFromOptions((*ctx)->rand)); 55 CHKERRQ(PetscDrawCreate(comm,host,label,x,y,m,n,&win)); 56 CHKERRQ(PetscDrawSetFromOptions(win)); 57 CHKERRQ(PetscDrawSPCreate(win,1,&(*ctx)->drawsp)); 58 CHKERRQ(KSPCreate(comm,&(*ctx)->ksp)); 59 CHKERRQ(KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_")); /* this is wrong, used use also prefix from the TS */ 60 CHKERRQ(KSPSetType((*ctx)->ksp,KSPGMRES)); 61 CHKERRQ(KSPGMRESSetRestart((*ctx)->ksp,200)); 62 CHKERRQ(KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200)); 63 CHKERRQ(KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE)); 64 CHKERRQ(KSPSetFromOptions((*ctx)->ksp)); 65 CHKERRQ(KSPGetPC((*ctx)->ksp,&pc)); 66 CHKERRQ(PCSetType(pc,PCNONE)); 67 68 (*ctx)->howoften = howoften; 69 (*ctx)->computeexplicitly = PETSC_FALSE; 70 71 CHKERRQ(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(0); 79 } 80 81 static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg) 82 { 83 PetscReal yr,yi; 84 85 PetscFunctionBegin; 86 CHKERRQ(TSComputeLinearStability(ts,xr,xi,&yr,&yi)); 87 if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE; 88 else *flg = PETSC_FALSE; 89 PetscFunctionReturn(0); 90 } 91 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(0); /* -1 indicates interpolated solution */ 105 if (!step) PetscFunctionReturn(0); 106 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 107 CHKERRQ(VecDuplicate(v,&xdot)); 108 CHKERRQ(TSGetSNES(ts,&snes)); 109 CHKERRQ(SNESGetJacobian(snes,&A,&B,NULL,NULL)); 110 CHKERRQ(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 CHKERRQ(SNESComputeJacobian(snes,v,A,B)); 119 120 ts->time_step = time_step_save; 121 122 CHKERRQ(KSPSetOperators(ksp,B,B)); 123 CHKERRQ(VecGetSize(v,&n)); 124 if (n < 200) its = n; 125 CHKERRQ(KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its)); 126 CHKERRQ(VecSetRandom(xdot,ctx->rand)); 127 CHKERRQ(KSPSolve(ksp,xdot,xdot)); 128 CHKERRQ(VecDestroy(&xdot)); 129 CHKERRQ(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 CHKERRQ(PetscDrawSPReset(drawsp)); 139 CHKERRQ(PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax)); 140 CHKERRQ(PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c)); 141 if (ctx->computeexplicitly) { 142 CHKERRQ(KSPComputeEigenvaluesExplicitly(ksp,n,r,c)); 143 neig = n; 144 } else { 145 CHKERRQ(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 CHKERRQ(TSComputeLinearStability(ts,r[i],c[i],&fr,&fi)); 153 if ((fr*fr + fi*fi) > 1.0) { 154 CHKERRQ(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))); 155 } 156 } 157 CHKERRQ(PetscDrawSPAddPoint(drawsp,r+i,c+i)); 158 } 159 CHKERRQ(PetscFree2(r,c)); 160 CHKERRQ(PetscDrawSPGetDraw(drawsp,&draw)); 161 CHKERRQ(PetscDrawGetPause(draw,&pause)); 162 CHKERRQ(PetscDrawSetPause(draw,0.0)); 163 CHKERRQ(PetscDrawSPDraw(drawsp,PETSC_TRUE)); 164 CHKERRQ(PetscDrawSetPause(draw,pause)); 165 if (ts->ops->linearstability) { 166 CHKERRQ(PetscDrawSPGetAxis(drawsp,&axis)); 167 CHKERRQ(PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax)); 168 CHKERRQ(PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts)); 169 CHKERRQ(PetscDrawSPDraw(drawsp,PETSC_FALSE)); 170 } 171 CHKERRQ(PetscDrawSPSave(drawsp)); 172 } 173 CHKERRQ(MatDestroy(&B)); 174 } 175 PetscFunctionReturn(0); 176 } 177 178 /*@C 179 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 180 181 Collective on TSMonitorSPEigCtx 182 183 Input Parameter: 184 . ctx - the monitor context 185 186 Level: intermediate 187 188 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 189 @*/ 190 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 191 { 192 PetscDraw draw; 193 194 PetscFunctionBegin; 195 CHKERRQ(PetscDrawSPGetDraw((*ctx)->drawsp,&draw)); 196 CHKERRQ(PetscDrawDestroy(&draw)); 197 CHKERRQ(PetscDrawSPDestroy(&(*ctx)->drawsp)); 198 CHKERRQ(KSPDestroy(&(*ctx)->ksp)); 199 CHKERRQ(PetscRandomDestroy(&(*ctx)->rand)); 200 CHKERRQ(PetscFree(*ctx)); 201 PetscFunctionReturn(0); 202 } 203