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