Actual source code: tseig.c

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

      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:   PetscReal   xmin, xmax, ymin, ymax;
     13: };

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

     18:   Collective

     20:   Input Parameters:
     21: + comm     - the communicator to share the monitor
     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        - the horizontal screen coordinates of the upper left coordinate of the window
     25: . y        - the vertical coordinates of the upper left coordinate of the window
     26: . m        - the screen width in pixels
     27: . n        - the screen height in pixels
     28: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

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

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

     36:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

    121:     ts->time_step = time_step_save;

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

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

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

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

    180:   Collective

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

    185:   Level: intermediate

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

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

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