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 = SNESComputeJacobian(snes,v,&A,&B,&structure);CHKERRQ(ierr); 91 /* ierr = TSComputeIJacobian(ts,ptime,v,xdot,1.0,&A,&B,&structure,PETSC_FALSE);CHKERRQ(ierr); */ 92 ierr = MatScale(A,-1.0);CHKERRQ(ierr); 93 if (A != B) { 94 ierr = MatScale(B,-1.0);CHKERRQ(ierr); 95 } 96 ierr = KSPSetOperators(ksp,A,B,structure);CHKERRQ(ierr); 97 ierr = VecSetRandom(xdot,ctx->rand);CHKERRQ(ierr); 98 ierr = KSPSolve(ksp,xdot,xdot);CHKERRQ(ierr); 99 ierr = VecDestroy(&xdot);CHKERRQ(ierr); 100 ierr = KSPGetIterationNumber(ksp,&nits);CHKERRQ(ierr); 101 n = nits+2; 102 103 if (nits) { 104 ierr = PetscMalloc2(n,PetscReal,&r,n,PetscReal,&c);CHKERRQ(ierr); 105 ierr = KSPComputeEigenvalues(ksp,n,r,c,&neig);CHKERRQ(ierr); 106 for (i=0; i<neig; i++) { 107 ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr); 108 } 109 ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr); 110 ierr = PetscFree2(r,c);CHKERRQ(ierr); 111 } 112 } 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "TSMonitorSPEigCtxDestroy" 118 /*@C 119 TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with TSMonitorSPEigCtxCreate(). 120 121 Collective on TSMonitorSPEigCtx 122 123 Input Parameter: 124 . ctx - the monitor context 125 126 Level: intermediate 127 128 .keywords: TS, monitor, line graph, destroy 129 130 .seealso: TSMonitorSPEigCtxCreate(), TSMonitorSet(), TSMonitorSPEig(); 131 @*/ 132 PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx) 133 { 134 PetscDraw draw; 135 PetscErrorCode ierr; 136 137 PetscFunctionBegin; 138 ierr = PetscDrawSPGetDraw((*ctx)->drawsp,&draw);CHKERRQ(ierr); 139 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 140 ierr = PetscDrawSPDestroy(&(*ctx)->drawsp);CHKERRQ(ierr); 141 ierr = KSPDestroy(&(*ctx)->ksp);CHKERRQ(ierr); 142 ierr = PetscRandomDestroy(&(*ctx)->rand);CHKERRQ(ierr); 143 ierr = PetscFree(*ctx);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147