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 (*ctx)->howoften = howoften; 73 (*ctx)->computeexplicitly = PETSC_FALSE; 74 ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,PETSC_NULL);CHKERRQ(ierr); 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 #undef __FUNCT__ 84 #define __FUNCT__ "TSLinearStabilityIndicator" 85 static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr,PetscReal xi,PetscBool *flg) 86 { 87 PetscErrorCode ierr; 88 PetscReal yr,yi; 89 90 PetscFunctionBegin; 91 ierr = TSComputeLinearStability(ts,xr,xi,&yr,&yi);CHKERRQ(ierr); 92 if ((yr*yr + yi*yi) <= 1.0) *flg = PETSC_TRUE; 93 else *flg = PETSC_FALSE; 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "TSMonitorSPEig" 99 PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 100 { 101 TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx; 102 PetscErrorCode ierr; 103 KSP ksp = ctx->ksp; 104 PetscInt n,N,nits,neig,i,its = 200; 105 PetscReal *r,*c,time_step_save; 106 PetscDrawSP drawsp = ctx->drawsp; 107 MatStructure structure; 108 Mat A,B; 109 Vec xdot; 110 SNES snes; 111 112 PetscFunctionBegin; 113 if (!step) PetscFunctionReturn(0); 114 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && (step == -1))){ 115 ierr = VecDuplicate(v,&xdot);CHKERRQ(ierr); 116 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 117 ierr = SNESGetJacobian(snes,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 118 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr); 119 /* 120 This doesn't work because methods keep and use internal information about the shift so it 121 seems we would need code for each method to trick the correct Jacobian in being computed. 122 */ 123 time_step_save = ts->time_step; 124 ts->time_step = PETSC_MAX_REAL; 125 ierr = SNESComputeJacobian(snes,v,&A,&B,&structure);CHKERRQ(ierr); 126 ts->time_step = time_step_save; 127 128 ierr = KSPSetOperators(ksp,B,B,structure);CHKERRQ(ierr); 129 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 130 if (n < 200) its = n; 131 ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);CHKERRQ(ierr); 132 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 133 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 134 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 135 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 136 N = nits+2; 137 138 if (nits) { 139 PetscDraw draw; 140 PetscReal pause; 141 PetscDrawAxis axis; 142 PetscReal xmin,xmax,ymin,ymax; 143 144 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr); 145 ierr = PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);CHKERRQ(ierr); 146 ierr = PetscMalloc2(PetscMax(n,N),PetscReal,&r,PetscMax(n,N),PetscReal,&c);CHKERRQ(ierr); 147 if (ctx->computeexplicitly) { 148 ierr = KSPComputeEigenvaluesExplicitly(ksp,n,r,c);CHKERRQ(ierr); 149 neig = n; 150 } else { 151 ierr = KSPComputeEigenvalues(ksp,N,r,c,&neig);CHKERRQ(ierr); 152 } 153 /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */ 154 for (i=0; i<neig; i++) r[i] = -r[i]; 155 for (i=0; i<neig; i++) { 156 if (ts->ops->linearstability) { 157 PetscReal fr,fi; 158 ierr = TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);CHKERRQ(ierr); 159 if ((fr*fr + fi*fi) > 1.0) { 160 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); 161 } 162 } 163 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 164 } 165 ierr = PetscFree2(r,c);CHKERRQ(ierr); 166 ierr = PetscDrawSPGetDraw(drawsp,&draw);CHKERRQ(ierr); 167 ierr = PetscDrawGetPause(draw,&pause);CHKERRQ(ierr); 168 ierr = PetscDrawSetPause(draw,0.0);CHKERRQ(ierr); 169 ierr = PetscDrawSPDraw(drawsp,PETSC_TRUE);CHKERRQ(ierr); 170 ierr = PetscDrawSetPause(draw,pause);CHKERRQ(ierr); 171 172 if (ts->ops->linearstability) { 173 ierr = PetscDrawSPGetAxis(drawsp,&axis);CHKERRQ(ierr); 174 ierr = PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);CHKERRQ(ierr); 175 ierr = PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);CHKERRQ(ierr); 176 ierr = PetscDrawSPDraw(drawsp,PETSC_FALSE);CHKERRQ(ierr); 177 } 178 } 179 ierr = MatDestroy(&B);CHKERRQ(ierr); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 186 /*@C 187 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 188 189 Collective on TSMonitorSPEigCtx 190 191 Input Parameter: 192 . ctx - the monitor context 193 194 Level: intermediate 195 196 .keywords: TS, monitor, line graph, destroy 197 198 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 199 @*/ 200 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 201 { 202 PetscDraw draw; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 207 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 208 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 209 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 210 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 211 ierr = PetscFree(*ctx);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214 215 216 217