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 PetscScalar shift; 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(struct _n_TSMonitorSPEigCtx,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 (*ctx)->howoften = howoften; 74 (*ctx)->computeexplicitly = PETSC_FALSE; 75 ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,PETSC_NULL);CHKERRQ(ierr); 76 (*ctx)->shift = 0.0; 77 ierr = PetscOptionsGetScalar(PETSC_NULL,"-ts_monitor_sp_eig_shift",&(*ctx)->shift,PETSC_NULL);CHKERRQ(ierr); 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; 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 ts->time_step = PETSC_MAX_REAL; 126 ierr = SNESComputeJacobian(snes,v,&A,&B,&structure);CHKERRQ(ierr); 127 */ 128 ierr = TSComputeIJacobian(ts,ptime,v,xdot,0.0,&B,&B,&structure,PETSC_FALSE);CHKERRQ(ierr); 129 ierr = MatScale(B,-1.0*ts->time_step);CHKERRQ(ierr); 130 ierr = MatShift(B,ctx->shift);CHKERRQ(ierr); 131 132 ierr = KSPSetOperators(ksp,B,B,structure);CHKERRQ(ierr); 133 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 134 if (n < 200) its = n; 135 ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);CHKERRQ(ierr); 136 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 137 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 138 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 139 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 140 N = nits+2; 141 142 if (nits) { 143 PetscDraw draw; 144 PetscReal pause; 145 PetscDrawAxis axis; 146 PetscReal xmin,xmax,ymin,ymax; 147 148 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr); 149 ierr = PetscDrawSPSetLimits(drawsp,ctx->xmin,ctx->xmax,ctx->ymin,ctx->ymax);CHKERRQ(ierr); 150 ierr = PetscMalloc2(PetscMax(n,N),PetscReal,&r,PetscMax(n,N),PetscReal,&c);CHKERRQ(ierr); 151 if (ctx->computeexplicitly) { 152 ierr = KSPComputeEigenvaluesExplicitly(ksp,n,r,c);CHKERRQ(ierr); 153 neig = n; 154 } else { 155 ierr = KSPComputeEigenvalues(ksp,N,r,c,&neig);CHKERRQ(ierr); 156 } 157 for (i=0; i<neig; i++) { 158 r[i] -= PetscAbsScalar(ctx->shift); 159 if (ts->ops->linearstability) { 160 PetscReal fr,fi; 161 ierr = TSComputeLinearStability(ts,r[i],c[i],&fr,&fi);CHKERRQ(ierr); 162 if ((fr*fr + fi*fi) > 1.0) { 163 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); 164 } 165 } 166 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 167 } 168 ierr = PetscFree2(r,c);CHKERRQ(ierr); 169 ierr = PetscDrawSPGetDraw(drawsp,&draw);CHKERRQ(ierr); 170 ierr = PetscDrawGetPause(draw,&pause);CHKERRQ(ierr); 171 ierr = PetscDrawSetPause(draw,0.0);CHKERRQ(ierr); 172 ierr = PetscDrawSPDraw(drawsp,PETSC_TRUE);CHKERRQ(ierr); 173 ierr = PetscDrawSetPause(draw,pause);CHKERRQ(ierr); 174 175 if (ts->ops->linearstability) { 176 ierr = PetscDrawSPGetAxis(drawsp,&axis);CHKERRQ(ierr); 177 ierr = PetscDrawAxisGetLimits(axis,&xmin,&xmax,&ymin,&ymax);CHKERRQ(ierr); 178 ierr = PetscDrawIndicatorFunction(draw,xmin,xmax,ymin,ymax,PETSC_DRAW_CYAN,(PetscErrorCode (*)(void*,PetscReal,PetscReal,PetscBool*))TSLinearStabilityIndicator,ts);CHKERRQ(ierr); 179 ierr = PetscDrawSPDraw(drawsp,PETSC_FALSE);CHKERRQ(ierr); 180 } 181 } 182 ierr = MatDestroy(&B);CHKERRQ(ierr); 183 } 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 189 /*@C 190 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 191 192 Collective on TSMonitorSPEigCtx 193 194 Input Parameter: 195 . ctx - the monitor context 196 197 Level: intermediate 198 199 .keywords: TS, monitor, line graph, destroy 200 201 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 202 @*/ 203 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 204 { 205 PetscDraw draw; 206 PetscErrorCode ierr; 207 208 PetscFunctionBegin; 209 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 210 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 211 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 212 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 213 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 214 ierr = PetscFree(*ctx);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 219 220