Actual source code: tseig.c


  2: #include <petsc/private/tsimpl.h>
  3: #include <petscdraw.h>

  5: /* ------------------------------------------------------------------------*/
  6: struct _n_TSMonitorSPEigCtx {
  7:   PetscDrawSP drawsp;
  8:   KSP         ksp;
  9:   PetscInt    howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
 10:   PetscBool   computeexplicitly;
 11:   MPI_Comm    comm;
 12:   PetscRandom rand;
 13:   PetscReal   xmin, xmax, ymin, ymax;
 14: };

 16: /*@C
 17:    TSMonitorSPEigCtxCreate - Creates a context for use with `TS` to monitor the eigenvalues of the linearized operator

 19:    Collective

 21:    Input Parameters:
 22: +  host - the X display to open, or null for the local machine
 23: .  label - the title to put in the title bar
 24: .  x, y - the screen coordinates of the upper left coordinate of the window
 25: .  m, n - the screen width and height in pixels
 26: -  howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

 28:    Output Parameter:
 29: .  ctx - the context

 31:    Options Database Key:
 32: .  -ts_monitor_sp_eig - plot egienvalues of linearized right hand side

 34:    Level: intermediate

 36:    Notes:
 37:    Use `TSMonitorSPEigCtxDestroy()` to destroy the context

 39:    Currently only works if the Jacobian is provided explicitly.

 41:    Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.

 43: .seealso: [](chapter_ts), `TSMonitorSPEigTimeStep()`, `TSMonitorSet()`, `TSMonitorLGSolution()`, `TSMonitorLGError()`
 44: @*/
 45: PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm comm, const char host[], const char label[], int x, int y, int m, int n, PetscInt howoften, TSMonitorSPEigCtx *ctx)
 46: {
 47:   PetscDraw win;
 48:   PC        pc;

 50:   PetscNew(ctx);
 51:   PetscRandomCreate(comm, &(*ctx)->rand);
 52:   PetscRandomSetFromOptions((*ctx)->rand);
 53:   PetscDrawCreate(comm, host, label, x, y, m, n, &win);
 54:   PetscDrawSetFromOptions(win);
 55:   PetscDrawSPCreate(win, 1, &(*ctx)->drawsp);
 56:   KSPCreate(comm, &(*ctx)->ksp);
 57:   KSPSetOptionsPrefix((*ctx)->ksp, "ts_monitor_sp_eig_"); /* this is wrong, used use also prefix from the TS */
 58:   KSPSetType((*ctx)->ksp, KSPGMRES);
 59:   KSPGMRESSetRestart((*ctx)->ksp, 200);
 60:   KSPSetTolerances((*ctx)->ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, 200);
 61:   KSPSetComputeSingularValues((*ctx)->ksp, PETSC_TRUE);
 62:   KSPSetFromOptions((*ctx)->ksp);
 63:   KSPGetPC((*ctx)->ksp, &pc);
 64:   PCSetType(pc, PCNONE);

 66:   (*ctx)->howoften          = howoften;
 67:   (*ctx)->computeexplicitly = PETSC_FALSE;

 69:   PetscOptionsGetBool(NULL, NULL, "-ts_monitor_sp_eig_explicitly", &(*ctx)->computeexplicitly, NULL);

 71:   (*ctx)->comm = comm;
 72:   (*ctx)->xmin = -2.1;
 73:   (*ctx)->xmax = 1.1;
 74:   (*ctx)->ymin = -1.1;
 75:   (*ctx)->ymax = 1.1;
 76:   return 0;
 77: }

 79: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr, PetscReal xi, PetscBool *flg)
 80: {
 81:   PetscReal yr, yi;

 83:   TSComputeLinearStability(ts, xr, xi, &yr, &yi);
 84:   if ((yr * yr + yi * yi) <= 1.0) *flg = PETSC_TRUE;
 85:   else *flg = PETSC_FALSE;
 86:   return 0;
 87: }

 89: PetscErrorCode TSMonitorSPEig(TS ts, PetscInt step, PetscReal ptime, Vec v, void *monctx)
 90: {
 91:   TSMonitorSPEigCtx ctx = (TSMonitorSPEigCtx)monctx;
 92:   KSP               ksp = ctx->ksp;
 93:   PetscInt          n, N, nits, neig, i, its = 200;
 94:   PetscReal        *r, *c, time_step_save;
 95:   PetscDrawSP       drawsp = ctx->drawsp;
 96:   Mat               A, B;
 97:   Vec               xdot;
 98:   SNES              snes;

100:   if (step < 0) return 0; /* -1 indicates interpolated solution */
101:   if (!step) return 0;
102:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
103:     VecDuplicate(v, &xdot);
104:     TSGetSNES(ts, &snes);
105:     SNESGetJacobian(snes, &A, &B, NULL, NULL);
106:     MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B);
107:     /*
108:        This doesn't work because methods keep and use internal information about the shift so it
109:        seems we would need code for each method to trick the correct Jacobian in being computed.
110:      */
111:     time_step_save = ts->time_step;
112:     ts->time_step  = PETSC_MAX_REAL;

114:     SNESComputeJacobian(snes, v, A, B);

116:     ts->time_step = time_step_save;

118:     KSPSetOperators(ksp, B, B);
119:     VecGetSize(v, &n);
120:     if (n < 200) its = n;
121:     KSPSetTolerances(ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, its);
122:     VecSetRandom(xdot, ctx->rand);
123:     KSPSolve(ksp, xdot, xdot);
124:     VecDestroy(&xdot);
125:     KSPGetIterationNumber(ksp, &nits);
126:     N = nits + 2;

128:     if (nits) {
129:       PetscDraw     draw;
130:       PetscReal     pause;
131:       PetscDrawAxis axis;
132:       PetscReal     xmin, xmax, ymin, ymax;

134:       PetscDrawSPReset(drawsp);
135:       PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax);
136:       PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c);
137:       if (ctx->computeexplicitly) {
138:         KSPComputeEigenvaluesExplicitly(ksp, n, r, c);
139:         neig = n;
140:       } else {
141:         KSPComputeEigenvalues(ksp, N, r, c, &neig);
142:       }
143:       /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
144:       for (i = 0; i < neig; i++) r[i] = -r[i];
145:       for (i = 0; i < neig; i++) {
146:         if (ts->ops->linearstability) {
147:           PetscReal fr, fi;
148:           TSComputeLinearStability(ts, r[i], c[i], &fr, &fi);
149:           if ((fr * fr + fi * fi) > 1.0) PetscPrintf(ctx->comm, "Linearized Eigenvalue %g + %g i linear stability function %g norm indicates unstable scheme \n", (double)r[i], (double)c[i], (double)(fr * fr + fi * fi));
150:         }
151:         PetscDrawSPAddPoint(drawsp, r + i, c + i);
152:       }
153:       PetscFree2(r, c);
154:       PetscDrawSPGetDraw(drawsp, &draw);
155:       PetscDrawGetPause(draw, &pause);
156:       PetscDrawSetPause(draw, 0.0);
157:       PetscDrawSPDraw(drawsp, PETSC_TRUE);
158:       PetscDrawSetPause(draw, pause);
159:       if (ts->ops->linearstability) {
160:         PetscDrawSPGetAxis(drawsp, &axis);
161:         PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax);
162:         PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode(*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts);
163:         PetscDrawSPDraw(drawsp, PETSC_FALSE);
164:       }
165:       PetscDrawSPSave(drawsp);
166:     }
167:     MatDestroy(&B);
168:   }
169:   return 0;
170: }

172: /*@C
173:    TSMonitorSPEigCtxDestroy - Destroys a scatter plot context that was created with `TSMonitorSPEigCtxCreate()`.

175:    Collective

177:    Input Parameter:
178: .  ctx - the monitor context

180:    Level: intermediate

182:    Note:
183:    Should be passed to `TSMonitorSet()` along with `TSMonitorSPEig()` an the context created with `TSMonitorSPEigCtxCreate()`

185: .seealso: [](chapter_ts), `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig();`
186: @*/
187: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
188: {
189:   PetscDraw draw;

191:   PetscDrawSPGetDraw((*ctx)->drawsp, &draw);
192:   PetscDrawDestroy(&draw);
193:   PetscDrawSPDestroy(&(*ctx)->drawsp);
194:   KSPDestroy(&(*ctx)->ksp);
195:   PetscRandomDestroy(&(*ctx)->rand);
196:   PetscFree(*ctx);
197:   return 0;
198: }