Actual source code: ex9adj.c

      1: static char help[] = "Basic equation for generator stability analysis.\n";

      3: /*F

      5: \begin{eqnarray}
      6:                  \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
      7:                  \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
      8: \end{eqnarray}

     10:   Ensemble of initial conditions
     11:    ./ex9 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly

     13:   Fault at .1 seconds
     14:    ./ex9 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly

     16:   Initial conditions same as when fault is ended
     17:    ./ex9 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly

     19: F*/

     21: /*
     22:    Include "petscts.h" so that we can use TS solvers.  Note that this
     23:    file automatically includes:
     24:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     25:      petscmat.h - matrices
     26:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     27:      petscviewer.h - viewers               petscpc.h  - preconditioners
     28:      petscksp.h   - linear solvers
     29: */

     31: #include <petscts.h>

     33: typedef struct {
     34:   PetscScalar H, D, omega_b, omega_s, Pmax, Pm, E, V, X, u_s, c;
     35:   PetscInt    beta;
     36:   PetscReal   tf, tcl;
     37: } AppCtx;

     39: PetscErrorCode PostStepFunction(TS ts)
     40: {
     41:   Vec                U;
     42:   PetscReal          t;
     43:   const PetscScalar *u;

     45:   PetscFunctionBegin;
     46:   PetscCall(TSGetTime(ts, &t));
     47:   PetscCall(TSGetSolution(ts, &U));
     48:   PetscCall(VecGetArrayRead(U, &u));
     49:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "delta(%3.2f) = %8.7f\n", (double)t, (double)u[0]));
     50:   PetscCall(VecRestoreArrayRead(U, &u));
     51:   PetscFunctionReturn(PETSC_SUCCESS);
     52: }

     54: /*
     55:      Defines the ODE passed to the ODE solver
     56: */
     57: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
     58: {
     59:   PetscScalar       *f, Pmax;
     60:   const PetscScalar *u;

     62:   PetscFunctionBegin;
     63:   /*  The next three lines allow us to access the entries of the vectors directly */
     64:   PetscCall(VecGetArrayRead(U, &u));
     65:   PetscCall(VecGetArray(F, &f));
     66:   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
     67:   else Pmax = ctx->Pmax;

     69:   f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
     70:   f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);

     72:   PetscCall(VecRestoreArrayRead(U, &u));
     73:   PetscCall(VecRestoreArray(F, &f));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: /*
     78:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
     79: */
     80: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
     81: {
     82:   PetscInt           rowcol[] = {0, 1};
     83:   PetscScalar        J[2][2], Pmax;
     84:   const PetscScalar *u;

     86:   PetscFunctionBegin;
     87:   PetscCall(VecGetArrayRead(U, &u));
     88:   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
     89:   else Pmax = ctx->Pmax;

     91:   J[0][0] = 0;
     92:   J[0][1] = ctx->omega_b;
     93:   J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
     94:   J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);

     96:   PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
     97:   PetscCall(VecRestoreArrayRead(U, &u));

     99:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    100:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    101:   if (A != B) {
    102:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    103:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    104:   }
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx0)
    109: {
    110:   PetscInt    row[] = {0, 1}, col[] = {0};
    111:   PetscScalar J[2][1];
    112:   AppCtx     *ctx = (AppCtx *)ctx0;

    114:   PetscFunctionBeginUser;
    115:   J[0][0] = 0;
    116:   J[1][0] = ctx->omega_s / (2.0 * ctx->H);
    117:   PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
    118:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    119:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
    124: {
    125:   PetscScalar       *r;
    126:   const PetscScalar *u;

    128:   PetscFunctionBegin;
    129:   PetscCall(VecGetArrayRead(U, &u));
    130:   PetscCall(VecGetArray(R, &r));
    131:   r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
    132:   PetscCall(VecRestoreArray(R, &r));
    133:   PetscCall(VecRestoreArrayRead(U, &u));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
    138: {
    139:   PetscScalar        ru[1];
    140:   const PetscScalar *u;
    141:   PetscInt           row[] = {0}, col[] = {0};

    143:   PetscFunctionBegin;
    144:   PetscCall(VecGetArrayRead(U, &u));
    145:   ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
    146:   PetscCall(VecRestoreArrayRead(U, &u));
    147:   PetscCall(MatSetValues(DRDU, 1, row, 1, col, ru, INSERT_VALUES));
    148:   PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
    149:   PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    151: }

    153: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, AppCtx *ctx)
    154: {
    155:   PetscFunctionBegin;
    156:   PetscCall(MatZeroEntries(DRDP));
    157:   PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
    158:   PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
    163: {
    164:   PetscScalar        sensip;
    165:   const PetscScalar *x, *y;

    167:   PetscFunctionBegin;
    168:   PetscCall(VecGetArrayRead(lambda, &x));
    169:   PetscCall(VecGetArrayRead(mu, &y));
    170:   sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
    171:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt parameter pm: %.7f \n", (double)sensip));
    172:   PetscCall(VecRestoreArrayRead(lambda, &x));
    173:   PetscCall(VecRestoreArrayRead(mu, &y));
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: int main(int argc, char **argv)
    178: {
    179:   TS           ts, quadts; /* ODE integrator */
    180:   Vec          U;          /* solution will be stored here */
    181:   Mat          A;          /* Jacobian matrix */
    182:   Mat          Jacp;       /* Jacobian matrix */
    183:   Mat          DRDU, DRDP;
    184:   PetscMPIInt  size;
    185:   PetscInt     n = 2;
    186:   AppCtx       ctx;
    187:   PetscScalar *u;
    188:   PetscReal    du[2]    = {0.0, 0.0};
    189:   PetscBool    ensemble = PETSC_FALSE, flg1, flg2;
    190:   PetscReal    ftime;
    191:   PetscInt     steps;
    192:   PetscScalar *x_ptr, *y_ptr;
    193:   Vec          lambda[1], q, mu[1];

    195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    196:      Initialize program
    197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    198:   PetscFunctionBeginUser;
    199:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    200:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    201:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    204:     Create necessary matrix and vectors
    205:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    206:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    207:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    208:   PetscCall(MatSetType(A, MATDENSE));
    209:   PetscCall(MatSetFromOptions(A));
    210:   PetscCall(MatSetUp(A));

    212:   PetscCall(MatCreateVecs(A, &U, NULL));

    214:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
    215:   PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
    216:   PetscCall(MatSetFromOptions(Jacp));
    217:   PetscCall(MatSetUp(Jacp));

    219:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &DRDP));
    220:   PetscCall(MatSetUp(DRDP));
    221:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, &DRDU));
    222:   PetscCall(MatSetUp(DRDU));

    224:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    225:     Set runtime options
    226:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    227:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
    228:   {
    229:     ctx.beta    = 2;
    230:     ctx.c       = 10000.0;
    231:     ctx.u_s     = 1.0;
    232:     ctx.omega_s = 1.0;
    233:     ctx.omega_b = 120.0 * PETSC_PI;
    234:     ctx.H       = 5.0;
    235:     PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
    236:     ctx.D = 5.0;
    237:     PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
    238:     ctx.E    = 1.1378;
    239:     ctx.V    = 1.0;
    240:     ctx.X    = 0.545;
    241:     ctx.Pmax = ctx.E * ctx.V / ctx.X;
    242:     PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
    243:     ctx.Pm = 1.1;
    244:     PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
    245:     ctx.tf  = 0.1;
    246:     ctx.tcl = 0.2;
    247:     PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
    248:     PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
    249:     PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
    250:     if (ensemble) {
    251:       ctx.tf  = -1;
    252:       ctx.tcl = -1;
    253:     }

    255:     PetscCall(VecGetArray(U, &u));
    256:     u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
    257:     u[1] = 1.0;
    258:     PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
    259:     n = 2;
    260:     PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
    261:     u[0] += du[0];
    262:     u[1] += du[1];
    263:     PetscCall(VecRestoreArray(U, &u));
    264:     if (flg1 || flg2) {
    265:       ctx.tf  = -1;
    266:       ctx.tcl = -1;
    267:     }
    268:   }
    269:   PetscOptionsEnd();

    271:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    272:      Create timestepping solver context
    273:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    274:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    275:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    276:   PetscCall(TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
    277:   PetscCall(TSSetType(ts, TSRK));
    278:   PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
    279:   PetscCall(TSSetRHSJacobian(ts, A, A, (TSRHSJacobianFn *)RHSJacobian, &ctx));
    280:   PetscCall(TSCreateQuadratureTS(ts, PETSC_TRUE, &quadts));
    281:   PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &ctx));
    282:   PetscCall(TSSetRHSJacobian(quadts, DRDU, DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, &ctx));
    283:   PetscCall(TSSetRHSJacobianP(quadts, DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, &ctx));
    284:   PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
    285:   PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &ctx));

    287:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    288:      Set initial conditions
    289:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    290:   PetscCall(TSSetSolution(ts, U));

    292:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    293:     Save trajectory of solution so that TSAdjointSolve() may be used
    294:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    295:   PetscCall(TSSetSaveTrajectory(ts));

    297:   PetscCall(MatCreateVecs(A, &lambda[0], NULL));
    298:   /*   Set initial conditions for the adjoint integration */
    299:   PetscCall(VecGetArray(lambda[0], &y_ptr));
    300:   y_ptr[0] = 0.0;
    301:   y_ptr[1] = 0.0;
    302:   PetscCall(VecRestoreArray(lambda[0], &y_ptr));

    304:   PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
    305:   PetscCall(VecGetArray(mu[0], &x_ptr));
    306:   x_ptr[0] = -1.0;
    307:   PetscCall(VecRestoreArray(mu[0], &x_ptr));

    309:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    310:      Set solver options
    311:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    312:   PetscCall(TSSetMaxTime(ts, 10.0));
    313:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    314:   PetscCall(TSSetTimeStep(ts, .01));
    315:   PetscCall(TSSetFromOptions(ts));

    317:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    318:      Solve nonlinear system
    319:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    320:   if (ensemble) {
    321:     for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
    322:       PetscCall(VecGetArray(U, &u));
    323:       u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
    324:       u[1] = ctx.omega_s;
    325:       u[0] += du[0];
    326:       u[1] += du[1];
    327:       PetscCall(VecRestoreArray(U, &u));
    328:       PetscCall(TSSetTimeStep(ts, .01));
    329:       PetscCall(TSSolve(ts, U));
    330:     }
    331:   } else {
    332:     PetscCall(TSSolve(ts, U));
    333:   }
    334:   PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
    335:   PetscCall(TSGetSolveTime(ts, &ftime));
    336:   PetscCall(TSGetStepNumber(ts, &steps));

    338:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    339:      Adjoint model starts here
    340:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    341:   /*   Set initial conditions for the adjoint integration */
    342:   PetscCall(VecGetArray(lambda[0], &y_ptr));
    343:   y_ptr[0] = 0.0;
    344:   y_ptr[1] = 0.0;
    345:   PetscCall(VecRestoreArray(lambda[0], &y_ptr));

    347:   PetscCall(VecGetArray(mu[0], &x_ptr));
    348:   x_ptr[0] = -1.0;
    349:   PetscCall(VecRestoreArray(mu[0], &x_ptr));

    351:   PetscCall(TSAdjointSolve(ts));

    353:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0]  d[Psi(tf)]/d[omega0]\n"));
    354:   PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
    355:   PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
    356:   PetscCall(TSGetCostIntegral(ts, &q));
    357:   PetscCall(VecView(q, PETSC_VIEWER_STDOUT_WORLD));
    358:   PetscCall(VecGetArray(q, &x_ptr));
    359:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
    360:   PetscCall(VecRestoreArray(q, &x_ptr));

    362:   PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));

    364:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    365:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    366:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    367:   PetscCall(MatDestroy(&A));
    368:   PetscCall(MatDestroy(&Jacp));
    369:   PetscCall(MatDestroy(&DRDU));
    370:   PetscCall(MatDestroy(&DRDP));
    371:   PetscCall(VecDestroy(&U));
    372:   PetscCall(VecDestroy(&lambda[0]));
    373:   PetscCall(VecDestroy(&mu[0]));
    374:   PetscCall(TSDestroy(&ts));
    375:   PetscCall(PetscFinalize());
    376:   return 0;
    377: }

    379: /*TEST

    381:    build:
    382:       requires: !complex

    384:    test:
    385:       args: -viewer_binary_skip_info -ts_adapt_type none

    387: TEST*/