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