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 PetscCall(PetscNew(ctx)); 53 PetscCall(PetscRandomCreate(comm,&(*ctx)->rand)); 54 PetscCall(PetscRandomSetFromOptions((*ctx)->rand)); 55 PetscCall(PetscDrawCreate(comm,host,label,x,y,m,n,&win)); 56 PetscCall(PetscDrawSetFromOptions(win)); 57 PetscCall(PetscDrawSPCreate(win,1,&(*ctx)->drawsp)); 58 PetscCall(KSPCreate(comm,&(*ctx)->ksp)); 59 PetscCall(KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_")); /* this is wrong, used use also prefix from the TS */ 60 PetscCall(KSPSetType((*ctx)->ksp,KSPGMRES)); 61 PetscCall(KSPGMRESSetRestart((*ctx)->ksp,200)); 62 PetscCall(KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200)); 63 PetscCall(KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE)); 64 PetscCall(KSPSetFromOptions((*ctx)->ksp)); 65 PetscCall(KSPGetPC((*ctx)->ksp,&pc)); 66 PetscCall(PCSetType(pc,PCNONE)); 67 68 (*ctx)->howoften = howoften; 69 (*ctx)->computeexplicitly = PETSC_FALSE; 70 71 PetscCall(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 PetscCall(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 PetscCall(VecDuplicate(v,&xdot)); 108 PetscCall(TSGetSNES(ts,&snes)); 109 PetscCall(SNESGetJacobian(snes,&A,&B,NULL,NULL)); 110 PetscCall(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 PetscCall(SNESComputeJacobian(snes,v,A,B)); 119 120 ts->time_step = time_step_save; 121 122 PetscCall(KSPSetOperators(ksp,B,B)); 123 PetscCall(VecGetSize(v,&n)); 124 if (n < 200) its = n; 125 PetscCall(KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its)); 126 PetscCall(VecSetRandom(xdot,ctx->rand)); 127 PetscCall(KSPSolve(ksp,xdot,xdot)); 128 PetscCall(VecDestroy(&xdot)); 129 PetscCall(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 PetscCall(PetscDrawSPReset(drawsp)); 139 PetscCall(PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax)); 140 PetscCall(PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c)); 141 if (ctx->computeexplicitly) { 142 PetscCall(KSPComputeEigenvaluesExplicitly(ksp,n,r,c)); 143 neig = n; 144 } else { 145 PetscCall(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 PetscCall(TSComputeLinearStability(ts,r[i],c[i],&fr,&fi)); 153 if ((fr*fr + fi*fi) > 1.0) { 154 PetscCall(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 PetscCall(PetscDrawSPAddPoint(drawsp,r+i,c+i)); 158 } 159 PetscCall(PetscFree2(r,c)); 160 PetscCall(PetscDrawSPGetDraw(drawsp,&draw)); 161 PetscCall(PetscDrawGetPause(draw,&pause)); 162 PetscCall(PetscDrawSetPause(draw,0.0)); 163 PetscCall(PetscDrawSPDraw(drawsp,PETSC_TRUE)); 164 PetscCall(PetscDrawSetPause(draw,pause)); 165 if (ts->ops->linearstability) { 166 PetscCall(PetscDrawSPGetAxis(drawsp,&axis)); 167 PetscCall(PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax)); 168 PetscCall(PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts)); 169 PetscCall(PetscDrawSPDraw(drawsp,PETSC_FALSE)); 170 } 171 PetscCall(PetscDrawSPSave(drawsp)); 172 } 173 PetscCall(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 PetscCall(PetscDrawSPGetDraw((*ctx)->drawsp,&draw)); 196 PetscCall(PetscDrawDestroy(&draw)); 197 PetscCall(PetscDrawSPDestroy(&(*ctx)->drawsp)); 198 PetscCall(KSPDestroy(&(*ctx)->ksp)); 199 PetscCall(PetscRandomDestroy(&(*ctx)->rand)); 200 PetscCall(PetscFree(*ctx)); 201 PetscFunctionReturn(0); 202 } 203