Actual source code: ex3sa.c

      1: static char help[] = "Adjoint and tangent linear sensitivity analysis of the 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: F*/

     12: /*
     13:   This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities.
     14:   It computes the sensitivities of an integral cost function
     15:   \int c*max(0,\theta(t)-u_s)^beta dt
     16:   w.r.t. initial conditions and the parameter P_m.
     17:   Backward Euler method is used for time integration.
     18:   The discontinuities are detected with TSEvent.
     19:  */

     21: #include <petscts.h>
     22: #include "ex3.h"

     24: int main(int argc, char **argv)
     25: {
     26:   TS           ts, quadts; /* ODE integrator */
     27:   Vec          U;          /* solution will be stored here */
     28:   PetscMPIInt  size;
     29:   PetscInt     n = 2;
     30:   AppCtx       ctx;
     31:   PetscScalar *u;
     32:   PetscReal    du[2]    = {0.0, 0.0};
     33:   PetscBool    ensemble = PETSC_FALSE, flg1, flg2;
     34:   PetscReal    ftime;
     35:   PetscInt     steps;
     36:   PetscScalar *x_ptr, *y_ptr, *s_ptr;
     37:   Vec          lambda[1], q, mu[1];
     38:   PetscInt     direction[2];
     39:   PetscBool    terminate[2];
     40:   Mat          qgrad;
     41:   Mat          sp; /* Forward sensitivity matrix */
     42:   SAMethod     sa;

     44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     45:      Initialize program
     46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     47:   PetscFunctionBeginUser;
     48:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     49:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     50:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

     52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     53:     Create necessary matrix and vectors
     54:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     55:   PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jac));
     56:   PetscCall(MatSetSizes(ctx.Jac, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
     57:   PetscCall(MatSetType(ctx.Jac, MATDENSE));
     58:   PetscCall(MatSetFromOptions(ctx.Jac));
     59:   PetscCall(MatSetUp(ctx.Jac));
     60:   PetscCall(MatCreateVecs(ctx.Jac, &U, NULL));
     61:   PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jacp));
     62:   PetscCall(MatSetSizes(ctx.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
     63:   PetscCall(MatSetFromOptions(ctx.Jacp));
     64:   PetscCall(MatSetUp(ctx.Jacp));
     65:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &ctx.DRDP));
     66:   PetscCall(MatSetUp(ctx.DRDP));
     67:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &ctx.DRDU));
     68:   PetscCall(MatSetUp(ctx.DRDU));

     70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     71:     Set runtime options
     72:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     73:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
     74:   {
     75:     ctx.beta    = 2;
     76:     ctx.c       = 10000.0;
     77:     ctx.u_s     = 1.0;
     78:     ctx.omega_s = 1.0;
     79:     ctx.omega_b = 120.0 * PETSC_PI;
     80:     ctx.H       = 5.0;
     81:     PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
     82:     ctx.D = 5.0;
     83:     PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
     84:     ctx.E        = 1.1378;
     85:     ctx.V        = 1.0;
     86:     ctx.X        = 0.545;
     87:     ctx.Pmax     = ctx.E * ctx.V / ctx.X;
     88:     ctx.Pmax_ini = ctx.Pmax;
     89:     PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
     90:     ctx.Pm = 1.1;
     91:     PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
     92:     ctx.tf  = 0.1;
     93:     ctx.tcl = 0.2;
     94:     PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
     95:     PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
     96:     PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
     97:     if (ensemble) {
     98:       ctx.tf  = -1;
     99:       ctx.tcl = -1;
    100:     }

    102:     PetscCall(VecGetArray(U, &u));
    103:     u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
    104:     u[1] = 1.0;
    105:     PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
    106:     n = 2;
    107:     PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
    108:     u[0] += du[0];
    109:     u[1] += du[1];
    110:     PetscCall(VecRestoreArray(U, &u));
    111:     if (flg1 || flg2) {
    112:       ctx.tf  = -1;
    113:       ctx.tcl = -1;
    114:     }
    115:     sa = SA_ADJ;
    116:     PetscCall(PetscOptionsEnum("-sa_method", "Sensitivity analysis method (adj or tlm)", "", SAMethods, (PetscEnum)sa, (PetscEnum *)&sa, NULL));
    117:   }
    118:   PetscOptionsEnd();

    120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121:      Create timestepping solver context
    122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    123:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    124:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    125:   PetscCall(TSSetType(ts, TSBEULER));
    126:   PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
    127:   PetscCall(TSSetRHSJacobian(ts, ctx.Jac, ctx.Jac, (TSRHSJacobianFn *)RHSJacobian, &ctx));

    129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    130:      Set initial conditions
    131:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    132:   PetscCall(TSSetSolution(ts, U));

    134:   /*   Set RHS JacobianP */
    135:   PetscCall(TSSetRHSJacobianP(ts, ctx.Jacp, RHSJacobianP, &ctx));

    137:   PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
    138:   PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &ctx));
    139:   PetscCall(TSSetRHSJacobian(quadts, ctx.DRDU, ctx.DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, &ctx));
    140:   PetscCall(TSSetRHSJacobianP(quadts, ctx.DRDP, DRDPJacobianTranspose, &ctx));
    141:   if (sa == SA_ADJ) {
    142:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    143:       Save trajectory of solution so that TSAdjointSolve() may be used
    144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    145:     PetscCall(TSSetSaveTrajectory(ts));
    146:     PetscCall(MatCreateVecs(ctx.Jac, &lambda[0], NULL));
    147:     PetscCall(MatCreateVecs(ctx.Jacp, &mu[0], NULL));
    148:     PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
    149:   }

    151:   if (sa == SA_TLM) {
    152:     PetscScalar val[2];
    153:     PetscInt    row[] = {0, 1}, col[] = {0};

    155:     PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &qgrad));
    156:     PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &sp));
    157:     PetscCall(TSForwardSetSensitivities(ts, 1, sp));
    158:     PetscCall(TSForwardSetSensitivities(quadts, 1, qgrad));
    159:     val[0] = 1. / PetscSqrtScalar(1. - (ctx.Pm / ctx.Pmax) * (ctx.Pm / ctx.Pmax)) / ctx.Pmax;
    160:     val[1] = 0.0;
    161:     PetscCall(MatSetValues(sp, 2, row, 1, col, val, INSERT_VALUES));
    162:     PetscCall(MatAssemblyBegin(sp, MAT_FINAL_ASSEMBLY));
    163:     PetscCall(MatAssemblyEnd(sp, MAT_FINAL_ASSEMBLY));
    164:   }

    166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    167:      Set solver options
    168:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    169:   PetscCall(TSSetMaxTime(ts, 1.0));
    170:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    171:   PetscCall(TSSetTimeStep(ts, 0.03125));
    172:   PetscCall(TSSetFromOptions(ts));

    174:   direction[0] = direction[1] = 1;
    175:   terminate[0] = terminate[1] = PETSC_FALSE;

    177:   PetscCall(TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&ctx));

    179:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    180:      Solve nonlinear system
    181:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    182:   if (ensemble) {
    183:     for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
    184:       PetscCall(VecGetArray(U, &u));
    185:       u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
    186:       u[1] = ctx.omega_s;
    187:       u[0] += du[0];
    188:       u[1] += du[1];
    189:       PetscCall(VecRestoreArray(U, &u));
    190:       PetscCall(TSSetTimeStep(ts, 0.03125));
    191:       PetscCall(TSSolve(ts, U));
    192:     }
    193:   } else {
    194:     PetscCall(TSSolve(ts, U));
    195:   }
    196:   PetscCall(TSGetSolveTime(ts, &ftime));
    197:   PetscCall(TSGetStepNumber(ts, &steps));

    199:   if (sa == SA_ADJ) {
    200:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    201:        Adjoint model starts here
    202:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    203:     /*   Set initial conditions for the adjoint integration */
    204:     PetscCall(VecGetArray(lambda[0], &y_ptr));
    205:     y_ptr[0] = 0.0;
    206:     y_ptr[1] = 0.0;
    207:     PetscCall(VecRestoreArray(lambda[0], &y_ptr));

    209:     PetscCall(VecGetArray(mu[0], &x_ptr));
    210:     x_ptr[0] = 0.0;
    211:     PetscCall(VecRestoreArray(mu[0], &x_ptr));

    213:     PetscCall(TSAdjointSolve(ts));

    215:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n lambda: d[Psi(tf)]/d[phi0]  d[Psi(tf)]/d[omega0]\n"));
    216:     PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
    217:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n mu: d[Psi(tf)]/d[pm]\n"));
    218:     PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
    219:     PetscCall(TSGetCostIntegral(ts, &q));
    220:     PetscCall(VecGetArray(q, &x_ptr));
    221:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
    222:     PetscCall(VecRestoreArray(q, &x_ptr));
    223:     PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
    224:     PetscCall(VecGetArray(mu[0], &x_ptr));
    225:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)x_ptr[0]));
    226:     PetscCall(VecRestoreArray(mu[0], &x_ptr));
    227:     PetscCall(VecDestroy(&lambda[0]));
    228:     PetscCall(VecDestroy(&mu[0]));
    229:   }
    230:   if (sa == SA_TLM) {
    231:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n trajectory sensitivity: d[phi(tf)]/d[pm]  d[omega(tf)]/d[pm]\n"));
    232:     PetscCall(MatView(sp, PETSC_VIEWER_STDOUT_WORLD));
    233:     PetscCall(TSGetCostIntegral(ts, &q));
    234:     PetscCall(VecGetArray(q, &s_ptr));
    235:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(s_ptr[0] - ctx.Pm)));
    236:     PetscCall(VecRestoreArray(q, &s_ptr));
    237:     PetscCall(MatDenseGetArray(qgrad, &s_ptr));
    238:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)s_ptr[0]));
    239:     PetscCall(MatDenseRestoreArray(qgrad, &s_ptr));
    240:     PetscCall(MatDestroy(&qgrad));
    241:     PetscCall(MatDestroy(&sp));
    242:   }
    243:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    244:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    245:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    246:   PetscCall(MatDestroy(&ctx.Jac));
    247:   PetscCall(MatDestroy(&ctx.Jacp));
    248:   PetscCall(MatDestroy(&ctx.DRDU));
    249:   PetscCall(MatDestroy(&ctx.DRDP));
    250:   PetscCall(VecDestroy(&U));
    251:   PetscCall(TSDestroy(&ts));
    252:   PetscCall(PetscFinalize());
    253:   return 0;
    254: }

    256: /*TEST

    258:    build:
    259:       requires: !complex !single

    261:    test:
    262:       args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu

    264:    test:
    265:       suffix: 2
    266:       args: -sa_method tlm -ts_type cn -pc_type lu

    268:    test:
    269:       suffix: 3
    270:       args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp

    272: TEST*/