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