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 PetscErrorCode ierr; 50 PC pc; 51 52 PetscFunctionBegin; 53 ierr = PetscNew(ctx);CHKERRQ(ierr); 54 ierr = PetscRandomCreate(comm,&(*ctx)->rand);CHKERRQ(ierr); 55 ierr = PetscRandomSetFromOptions((*ctx)->rand);CHKERRQ(ierr); 56 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 57 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 58 ierr = PetscDrawSPCreate(win,1,&(*ctx)->drawsp);CHKERRQ(ierr); 59 ierr = KSPCreate(comm,&(*ctx)->ksp);CHKERRQ(ierr); 60 ierr = KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_");CHKERRQ(ierr); /* this is wrong, used use also prefix from the TS */ 61 ierr = KSPSetType((*ctx)->ksp,KSPGMRES);CHKERRQ(ierr); 62 ierr = KSPGMRESSetRestart((*ctx)->ksp,200);CHKERRQ(ierr); 63 ierr = KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);CHKERRQ(ierr); 64 ierr = KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);CHKERRQ(ierr); 65 ierr = KSPSetFromOptions((*ctx)->ksp);CHKERRQ(ierr); 66 ierr = KSPGetPC((*ctx)->ksp,&pc);CHKERRQ(ierr); 67 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 68 69 (*ctx)->howoften = howoften; 70 (*ctx)->computeexplicitly = PETSC_FALSE; 71 72 ierr = PetscOptionsGetBool(NULL,NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);CHKERRQ(ierr); 73 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 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg) 83 { 84 PetscErrorCode ierr; 85 PetscReal yr,yi; 86 87 PetscFunctionBegin; 88 ierr = TSComputeLinearStability(ts,xr,xi,&yr,&yi);CHKERRQ(ierr); 89 if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE; 90 else *flg = PETSC_FALSE; 91 PetscFunctionReturn(0); 92 } 93 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; 105 106 PetscFunctionBegin; 107 if (step < 0) PetscFunctionReturn(0); /* -1 indicates interpolated solution */ 108 if (!step) PetscFunctionReturn(0); 109 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 110 ierr = VecDuplicate(v,&xdot);CHKERRQ(ierr); 111 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 112 ierr = SNESGetJacobian(snes,&A,&B,NULL,NULL);CHKERRQ(ierr); 113 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr); 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; 120 121 ierr = SNESComputeJacobian(snes,v,A,B);CHKERRQ(ierr); 122 123 ts->time_step = time_step_save; 124 125 ierr = KSPSetOperators(ksp,B,B);CHKERRQ(ierr); 126 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 127 if (n < 200) its = n; 128 ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);CHKERRQ(ierr); 129 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 130 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 131 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 132 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 133 N = nits+2; 134 135 if (nits) { 136 PetscDraw draw; 137 PetscReal pause; 138 PetscDrawAxis axis; 139 PetscReal xmin,xmax,ymin,ymax; 140 141 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr); 142 ierr = PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);CHKERRQ(ierr); 143 ierr = PetscMalloc2(PetscMax(n,N),&r,PetscMax(n,N),&c);CHKERRQ(ierr); 144 if (ctx->computeexplicitly) { 145 ierr = KSPComputeEigenvaluesExplicitly(ksp,n,r,c);CHKERRQ(ierr); 146 neig = n; 147 } else { 148 ierr = KSPComputeEigenvalues(ksp,N,r,c,&neig);CHKERRQ(ierr); 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 ierr = TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);CHKERRQ(ierr); 156 if ((fr*fr + fi*fi) > 1.0) { 157 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); 158 } 159 } 160 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 161 } 162 ierr = PetscFree2(r,c);CHKERRQ(ierr); 163 ierr = PetscDrawSPGetDraw(drawsp,&draw);CHKERRQ(ierr); 164 ierr = PetscDrawGetPause(draw,&pause);CHKERRQ(ierr); 165 ierr = PetscDrawSetPause(draw,0.0);CHKERRQ(ierr); 166 ierr = PetscDrawSPDraw(drawsp,PETSC_TRUE);CHKERRQ(ierr); 167 ierr = PetscDrawSetPause(draw,pause);CHKERRQ(ierr); 168 if (ts->ops->linearstability) { 169 ierr = PetscDrawSPGetAxis(drawsp,&axis);CHKERRQ(ierr); 170 ierr = PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);CHKERRQ(ierr); 171 ierr = PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);CHKERRQ(ierr); 172 ierr = PetscDrawSPDraw(drawsp,PETSC_FALSE);CHKERRQ(ierr); 173 } 174 ierr = PetscDrawSPSave(drawsp);CHKERRQ(ierr); 175 } 176 ierr = MatDestroy(&B);CHKERRQ(ierr); 177 } 178 PetscFunctionReturn(0); 179 } 180 181 /*@C 182 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 183 184 Collective on TSMonitorSPEigCtx 185 186 Input Parameter: 187 . ctx - the monitor context 188 189 Level: intermediate 190 191 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 192 @*/ 193 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 194 { 195 PetscDraw draw; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 200 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 201 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 202 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 203 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 204 ierr = PetscFree(*ctx);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208