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