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 /* 103 This doesn't work because methods keep and use internal information about the shift so it 104 seems we would need code for each method to trick the correct Jacobian in being computed. 105 ts->time_step = PETSC_MAX_REAL; 106 ierr = SNESComputeJacobian(snes,v,&A,&B,&structure);CHKERRQ(ierr); 107 */ 108 ierr = TSComputeIJacobian(ts,ptime,v,xdot,0.0,&A,&B,&structure,PETSC_FALSE);CHKERRQ(ierr); 109 ierr = MatScale(A,-1.0);CHKERRQ(ierr); 110 ierr = MatShift(A,ctx->shift);CHKERRQ(ierr); 111 ierr = KSPSetOperators(ksp,A,A,structure);CHKERRQ(ierr); 112 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 113 if (n < 200) its = n; 114 ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,its);CHKERRQ(ierr); 115 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 116 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 117 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 118 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 119 N = nits+2; 120 121 if (nits) { 122 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr); 123 ierr = PetscMalloc2(N,PetscReal,&r,N,PetscReal,&c);CHKERRQ(ierr); 124 if (ctx->computeexplicitly) { 125 ierr = KSPComputeEigenvaluesExplicitly(ksp,n,r,c);CHKERRQ(ierr); 126 neig = n; 127 } else { 128 ierr = KSPComputeEigenvalues(ksp,N,r,c,&neig);CHKERRQ(ierr); 129 } 130 for (i=0; i<neig; i++) { 131 r[i] -= ctx->shift; 132 ierr = PetscPrintf(ctx->comm,"%g + %g i\n",r[i],c[i]);CHKERRQ(ierr); 133 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 134 } 135 ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr); 136 ierr = PetscFree2(r,c);CHKERRQ(ierr); 137 } 138 ierr = MatShift(A,-ctx->shift);CHKERRQ(ierr); 139 ierr = MatScale(A,-1.0);CHKERRQ(ierr); 140 } 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 146 /*@C 147 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 148 149 Collective on TSMonitorSPEigCtx 150 151 Input Parameter: 152 . ctx - the monitor context 153 154 Level: intermediate 155 156 .keywords: TS, monitor, line graph, destroy 157 158 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 159 @*/ 160 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 161 { 162 PetscDraw draw; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 167 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 168 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 169 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 170 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 171 ierr = PetscFree(*ctx);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175