1 2 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 3 4 /* ------------------------------------------------------------------------*/ 5 struct _n_TSMonitorSPEigCtx { 6 PetscDrawSP drawsp; 7 KSP ksp; 8 PetscRandom rand; 9 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 10 }; 11 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "TSMonitorSPEigCtxCreate" 15 /*@C 16 TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator 17 18 Collective on TS 19 20 Input Parameters: 21 + host - the X display to open, or null for the local machine 22 . label - the title to put in the title bar 23 . x, y - the screen coordinates of the upper left coordinate of the window 24 . m, n - the screen width and height in pixels 25 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 26 27 Output Parameter: 28 . ctx - the context 29 30 Options Database Key: 31 . -ts_monitor_eig - plot egienvalues of linearized right hand side 32 33 Notes: 34 Use TSMonitorSPEigCtxDestroy() to destroy. 35 36 Level: intermediate 37 38 .keywords: TS, monitor, line graph, residual, seealso 39 40 .seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 41 42 @*/ 43 PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx *ctx) 44 { 45 PetscDraw win; 46 PetscErrorCode ierr; 47 PC pc; 48 49 PetscFunctionBegin; 50 ierr = PetscNew(struct _n_TSMonitorSPEigCtx,ctx);CHKERRQ(ierr); 51 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 52 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 53 ierr = PetscDrawSPCreate(win,1,&(*ctx)->drawsp);CHKERRQ(ierr); 54 ierr = PetscRandomCreate(comm,&(*ctx)->rand);CHKERRQ(ierr); 55 ierr = PetscRandomSetFromOptions((*ctx)->rand);CHKERRQ(ierr); 56 ierr = KSPCreate(comm,&(*ctx)->ksp);CHKERRQ(ierr); 57 ierr = KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_eig_");CHKERRQ(ierr); /* this is wrong, used use also prefix from the TS */ 58 ierr = KSPSetType((*ctx)->ksp,KSPGMRES);CHKERRQ(ierr); 59 ierr = KSPGMRESSetRestart((*ctx)->ksp,200);CHKERRQ(ierr); 60 ierr = KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);CHKERRQ(ierr); 61 ierr = KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);CHKERRQ(ierr); 62 ierr = KSPSetFromOptions((*ctx)->ksp);CHKERRQ(ierr); 63 ierr = KSPGetPC((*ctx)->ksp,&pc);CHKERRQ(ierr); 64 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 65 (*ctx)->howoften = howoften; 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "TSMonitorSPEig" 71 PetscErrorCode TSMonitorSPEig(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 72 { 73 TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx) monctx; 74 PetscErrorCode ierr; 75 KSP ksp = ctx->ksp; 76 PetscInt n,nits,neig,i; 77 PetscReal *r,*c; 78 PetscDrawSP drawsp = ctx->drawsp; 79 MatStructure structure; 80 Mat A,B; 81 Vec xdot; 82 SNES snes; 83 84 PetscFunctionBegin; 85 if (!step) PetscFunctionReturn(0); 86 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && (step == -1))){ 87 ierr = VecDuplicate(v,&xdot);CHKERRQ(ierr); 88 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 89 ierr = SNESGetJacobian(snes,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 90 ierr = TSComputeIJacobian(ts,ptime,v,xdot,1.0,&A,&B,&structure,PETSC_FALSE);CHKERRQ(ierr); 91 ierr = KSPSetOperators(ksp,A,B,structure);CHKERRQ(ierr); 92 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 93 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 94 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 95 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 96 n = nits+2; 97 98 if (nits) { 99 ierr = PetscMalloc2(n,PetscReal,&r,n,PetscReal,&c);CHKERRQ(ierr); 100 ierr = KSPComputeEigenvalues(ksp,n,r,c,&neig);CHKERRQ(ierr); 101 for (i=0; i<neig; i++) { 102 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 103 } 104 ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr); 105 ierr = PetscFree2(r,c);CHKERRQ(ierr); 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 113 /*@C 114 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 115 116 Collective on TSMonitorSPEigCtx 117 118 Input Parameter: 119 . ctx - the monitor context 120 121 Level: intermediate 122 123 .keywords: TS, monitor, line graph, destroy 124 125 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 126 @*/ 127 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 128 { 129 PetscDraw draw; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 134 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 135 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 136 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 137 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 138 ierr = PetscFree(*ctx);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142