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: +  comm - the communicator to share the monitor
 23: .  host - the X display to open, or `NULL` for the local machine
 24: .  label - the title to put in the title bar
 25: .  x - the horizontal screen coordinates of the upper left coordinate of the window
 26: .  y - the vertical coordinates of the upper left coordinate of the window
 27: .  m - the screen width in pixels
 28: .  n - the screen height in pixels
 29: -  howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

 31:    Output Parameter:
 32: .  ctx - the context

 34:    Options Database Key:
 35: .  -ts_monitor_sp_eig - plot egienvalues of linearized right hand side

 37:    Level: intermediate

 39:    Notes:
 40:    Use `TSMonitorSPEigCtxDestroy()` to destroy the context

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

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

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

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

 70:   (*ctx)->howoften          = howoften;
 71:   (*ctx)->computeexplicitly = PETSC_FALSE;

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

 75:   (*ctx)->comm = comm;
 76:   (*ctx)->xmin = -2.1;
 77:   (*ctx)->xmax = 1.1;
 78:   (*ctx)->ymin = -1.1;
 79:   (*ctx)->ymax = 1.1;
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: static PetscErrorCode TSLinearStabilityIndicator(TS ts, PetscReal xr, PetscReal xi, PetscBool *flg)
 84: {
 85:   PetscReal yr, yi;

 87:   PetscFunctionBegin;
 88:   PetscCall(TSComputeLinearStability(ts, xr, xi, &yr, &yi));
 89:   if ((yr * yr + yi * yi) <= 1.0) *flg = PETSC_TRUE;
 90:   else *flg = PETSC_FALSE;
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

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

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

120:     PetscCall(SNESComputeJacobian(snes, v, A, B));

122:     ts->time_step = time_step_save;

124:     PetscCall(KSPSetOperators(ksp, B, B));
125:     PetscCall(VecGetSize(v, &n));
126:     if (n < 200) its = n;
127:     PetscCall(KSPSetTolerances(ksp, 1.e-10, PETSC_DEFAULT, PETSC_DEFAULT, its));
128:     PetscCall(VecSetRandom(xdot, ctx->rand));
129:     PetscCall(KSPSolve(ksp, xdot, xdot));
130:     PetscCall(VecDestroy(&xdot));
131:     PetscCall(KSPGetIterationNumber(ksp, &nits));
132:     N = nits + 2;

134:     if (nits) {
135:       PetscDraw     draw;
136:       PetscReal     pause;
137:       PetscDrawAxis axis;
138:       PetscReal     xmin, xmax, ymin, ymax;

140:       PetscCall(PetscDrawSPReset(drawsp));
141:       PetscCall(PetscDrawSPSetLimits(drawsp, ctx->xmin, ctx->xmax, ctx->ymin, ctx->ymax));
142:       PetscCall(PetscMalloc2(PetscMax(n, N), &r, PetscMax(n, N), &c));
143:       if (ctx->computeexplicitly) {
144:         PetscCall(KSPComputeEigenvaluesExplicitly(ksp, n, r, c));
145:         neig = n;
146:       } else {
147:         PetscCall(KSPComputeEigenvalues(ksp, N, r, c, &neig));
148:       }
149:       /* We used the positive operator to be able to reuse KSPs that require positive definiteness, now flip the spectrum as is conventional for ODEs */
150:       for (i = 0; i < neig; i++) r[i] = -r[i];
151:       for (i = 0; i < neig; i++) {
152:         if (ts->ops->linearstability) {
153:           PetscReal fr, fi;
154:           PetscCall(TSComputeLinearStability(ts, r[i], c[i], &fr, &fi));
155:           if ((fr * fr + fi * fi) > 1.0) PetscCall(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)));
156:         }
157:         PetscCall(PetscDrawSPAddPoint(drawsp, r + i, c + i));
158:       }
159:       PetscCall(PetscFree2(r, c));
160:       PetscCall(PetscDrawSPGetDraw(drawsp, &draw));
161:       PetscCall(PetscDrawGetPause(draw, &pause));
162:       PetscCall(PetscDrawSetPause(draw, 0.0));
163:       PetscCall(PetscDrawSPDraw(drawsp, PETSC_TRUE));
164:       PetscCall(PetscDrawSetPause(draw, pause));
165:       if (ts->ops->linearstability) {
166:         PetscCall(PetscDrawSPGetAxis(drawsp, &axis));
167:         PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, &ymin, &ymax));
168:         PetscCall(PetscDrawIndicatorFunction(draw, xmin, xmax, ymin, ymax, PETSC_DRAW_CYAN, (PetscErrorCode(*)(void *, PetscReal, PetscReal, PetscBool *))TSLinearStabilityIndicator, ts));
169:         PetscCall(PetscDrawSPDraw(drawsp, PETSC_FALSE));
170:       }
171:       PetscCall(PetscDrawSPSave(drawsp));
172:     }
173:     PetscCall(MatDestroy(&B));
174:   }
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

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

181:    Collective

183:    Input Parameter:
184: .  ctx - the monitor context

186:    Level: intermediate

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

191: .seealso: [](ch_ts), `TSMonitorSPEigCtxCreate()`, `TSMonitorSet()`, `TSMonitorSPEig();`
192: @*/
193: PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx *ctx)
194: {
195:   PetscDraw draw;

197:   PetscFunctionBegin;
198:   PetscCall(PetscDrawSPGetDraw((*ctx)->drawsp, &draw));
199:   PetscCall(PetscDrawDestroy(&draw));
200:   PetscCall(PetscDrawSPDestroy(&(*ctx)->drawsp));
201:   PetscCall(KSPDestroy(&(*ctx)->ksp));
202:   PetscCall(PetscRandomDestroy(&(*ctx)->rand));
203:   PetscCall(PetscFree(*ctx));
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }