1 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 }; 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)->shift = 0.0; 76 ierr = PetscOptionsGetScalar(PETSC_NULL,"-ts_monitor_sp_eig_shift",&(*ctx)->shift,PETSC_NULL);CHKERRQ(ierr); 77 (*ctx)->comm = comm; 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "TSMonitorSPEig" 83 PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 84 { 85 TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx; 86 PetscErrorCode ierr; 87 KSP ksp = ctx->ksp; 88 PetscInt n,N,nits,neig,i,its = 200; 89 PetscReal *r,*c; 90 PetscDrawSP drawsp = ctx->drawsp; 91 MatStructure structure; 92 Mat A,B; 93 Vec xdot; 94 SNES snes; 95 96 PetscFunctionBegin; 97 if (!step) PetscFunctionReturn(0); 98 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && (step == -1))){ 99 ierr = VecDuplicate(v,&xdot);CHKERRQ(ierr); 100 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 101 ierr = SNESGetJacobian(snes,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 102 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr); 103 /* 104 This doesn't work because methods keep and use internal information about the shift so it 105 seems we would need code for each method to trick the correct Jacobian in being computed. 106 ts->time_step = PETSC_MAX_REAL; 107 ierr = SNESComputeJacobian(snes,v,&A,&B,&structure);CHKERRQ(ierr); 108 */ 109 ierr = TSComputeIJacobian(ts,ptime,v,xdot,0.0,&B,&B,&structure,PETSC_FALSE);CHKERRQ(ierr); 110 ierr = MatScale(B,-1.0*ts->time_step);CHKERRQ(ierr); 111 ierr = MatShift(B,ctx->shift);CHKERRQ(ierr); 112 113 ierr = KSPSetOperators(ksp,B,B,structure);CHKERRQ(ierr); 114 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 115 if (n < 200) its = n; 116 ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);CHKERRQ(ierr); 117 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 118 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 119 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 120 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 121 N = nits+2; 122 123 if (nits) { 124 PetscDraw draw; 125 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr); 126 ierr = PetscDrawSPSetLimits(drawsp,-2.1,.1,-1.1,1.1);CHKERRQ(ierr); 127 ierr = PetscMalloc2(PetscMax(n,N),PetscReal,&r,PetscMax(n,N),PetscReal,&c);CHKERRQ(ierr); 128 if (ctx->computeexplicitly) { 129 ierr = KSPComputeEigenvaluesExplicitly(ksp,n,r,c);CHKERRQ(ierr); 130 neig = n; 131 } else { 132 ierr = KSPComputeEigenvalues(ksp,N,r,c,&neig);CHKERRQ(ierr); 133 } 134 for (i=0; i<neig; i++) { 135 r[i] -= ctx->shift; 136 ierr = PetscPrintf(ctx->comm,"%g + %g i\n",r[i],c[i]);CHKERRQ(ierr); 137 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 138 } 139 ierr = PetscDrawSPDraw(drawsp,PETSC_TRUE);CHKERRQ(ierr); 140 ierr = PetscDrawSPGetDraw(drawsp,&draw);CHKERRQ(ierr); 141 ierr = PetscDrawEllipse(draw,-1.0,0.0,2.0,2.0,PETSC_DRAW_CYAN);CHKERRQ(ierr); 142 ierr = PetscDrawSPDraw(drawsp,PETSC_FALSE);CHKERRQ(ierr); 143 ierr = PetscFree2(r,c);CHKERRQ(ierr); 144 } 145 ierr = MatDestroy(&B);CHKERRQ(ierr); 146 } 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 152 /*@C 153 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 154 155 Collective on TSMonitorSPEigCtx 156 157 Input Parameter: 158 . ctx - the monitor context 159 160 Level: intermediate 161 162 .keywords: TS, monitor, line graph, destroy 163 164 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 165 @*/ 166 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 167 { 168 PetscDraw draw; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 173 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 174 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 175 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 176 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 177 ierr = PetscFree(*ctx);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181