Actual source code: ex1fwd.c

      1: static char help[] = "Trajectory sensitivity of a hybrid system with state-dependent switchings.\n";

      3: /*
      4:   The dynamics is described by the ODE
      5:                   u_t = A_i u

      7:   where A_1 = [ 1  -100
      8:                 10  1  ],
      9:         A_2 = [ 1    10
     10:                -100  1 ].
     11:   The index i changes from 1 to 2 when u[1]=2.75u[0] and from 2 to 1 when u[1]=0.36u[0].
     12:   Initially u=[0 1]^T and i=1.

     14:   References:
     15: + * - H. Zhang, S. Abhyankar, E. Constantinescu, M. Mihai, Discrete Adjoint Sensitivity Analysis of Hybrid Dynamical Systems With Switching,
     16:       IEEE Transactions on Circuits and Systems I: Regular Papers, 64(5), May 2017
     17: - * - I. A. Hiskens, M.A. Pai, Trajectory Sensitivity Analysis of Hybrid Systems, IEEE Transactions on Circuits and Systems, Vol 47, No 2, February 2000
     18: */

     20: #include <petscts.h>

     22: typedef struct {
     23:   PetscScalar lambda1;
     24:   PetscScalar lambda2;
     25:   PetscInt    mode; /* mode flag*/
     26:   PetscReal   print_time;
     27: } AppCtx;

     29: PetscErrorCode MyMonitor(TS ts, PetscInt stepnum, PetscReal time, Vec U, void *ctx)
     30: {
     31:   AppCtx      *actx = (AppCtx *)ctx;
     32:   Mat          sp;
     33:   PetscScalar *u;
     34:   PetscInt     nump;
     35:   FILE        *f;

     37:   PetscFunctionBegin;
     38:   if (time >= actx->print_time) {
     39:     actx->print_time += 1. / 256.;
     40:     PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
     41:     PetscCall(MatDenseGetColumn(sp, 2, &u));
     42:     f = fopen("fwd_sp.out", "a");
     43:     PetscCall(PetscFPrintf(PETSC_COMM_WORLD, f, "%20.15lf %20.15lf %20.15lf\n", (double)time, (double)PetscRealPart(u[0]), (double)PetscRealPart(u[1])));
     44:     PetscCall(MatDenseRestoreColumn(sp, &u));
     45:     fclose(f);
     46:   }
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

     50: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx)
     51: {
     52:   AppCtx            *actx = (AppCtx *)ctx;
     53:   const PetscScalar *u;

     55:   PetscFunctionBegin;
     56:   PetscCall(VecGetArrayRead(U, &u));
     57:   if (actx->mode == 1) {
     58:     fvalue[0] = PetscRealPart(u[1] - actx->lambda1 * u[0]);
     59:   } else if (actx->mode == 2) {
     60:     fvalue[0] = PetscRealPart(u[1] - actx->lambda2 * u[0]);
     61:   }
     62:   PetscCall(VecRestoreArrayRead(U, &u));
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: PetscErrorCode ShiftGradients(TS ts, Vec U, AppCtx *actx)
     67: {
     68:   Mat          sp;
     69:   PetscScalar *x;
     70:   PetscScalar *u;
     71:   PetscScalar  tmp[2], A1[2][2], A2[2], denorm;
     72:   PetscInt     nump;

     74:   PetscFunctionBegin;
     75:   PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
     76:   PetscCall(VecGetArray(U, &u));

     78:   if (actx->mode == 1) {
     79:     denorm   = -actx->lambda1 * (u[0] - 100. * u[1]) + 1. * (10. * u[0] + u[1]);
     80:     A1[0][0] = 110. * u[1] * (-actx->lambda1) / denorm + 1.;
     81:     A1[1][0] = -110. * u[0] * (-actx->lambda1) / denorm;
     82:     A1[0][1] = 110. * u[1] * 1. / denorm;
     83:     A1[1][1] = -110. * u[0] * 1. / denorm + 1.;

     85:     A2[0] = 110. * u[1] * (-u[0]) / denorm;
     86:     A2[1] = -110. * u[0] * (-u[0]) / denorm;
     87:   } else {
     88:     denorm   = -actx->lambda2 * (u[0] + 10. * u[1]) + 1. * (-100. * u[0] + u[1]);
     89:     A1[0][0] = 110. * u[1] * (actx->lambda2) / denorm + 1;
     90:     A1[1][0] = -110. * u[0] * (actx->lambda2) / denorm;
     91:     A1[0][1] = -110. * u[1] * 1. / denorm;
     92:     A1[1][1] = 110. * u[0] * 1. / denorm + 1.;

     94:     A2[0] = 0;
     95:     A2[1] = 0;
     96:   }

     98:   PetscCall(VecRestoreArray(U, &u));

    100:   PetscCall(MatDenseGetColumn(sp, 0, &x));
    101:   tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
    102:   tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
    103:   x[0]   = tmp[0];
    104:   x[1]   = tmp[1];
    105:   PetscCall(MatDenseRestoreColumn(sp, &x));

    107:   PetscCall(MatDenseGetColumn(sp, 1, &x));
    108:   tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
    109:   tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
    110:   x[0]   = tmp[0];
    111:   x[1]   = tmp[1];
    112:   PetscCall(MatDenseRestoreColumn(sp, &x));

    114:   PetscCall(MatDenseGetColumn(sp, 2, &x));
    115:   tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
    116:   tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
    117:   x[0]   = tmp[0] + A2[0];
    118:   x[1]   = tmp[1] + A2[1];
    119:   PetscCall(MatDenseRestoreColumn(sp, &x));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
    124: {
    125:   AppCtx *actx = (AppCtx *)ctx;

    127:   PetscFunctionBegin;
    128:   /* PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); */
    129:   PetscCall(ShiftGradients(ts, U, actx));
    130:   if (actx->mode == 1) {
    131:     actx->mode = 2;
    132:     /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 1 to 2 at t = %f \n",(double)t)); */
    133:   } else if (actx->mode == 2) {
    134:     actx->mode = 1;
    135:     /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 2 to 1 at t = %f \n",(double)t)); */
    136:   }
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: /*
    141:      Defines the ODE passed to the ODE solver
    142: */
    143: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
    144: {
    145:   AppCtx            *actx = (AppCtx *)ctx;
    146:   PetscScalar       *f;
    147:   const PetscScalar *u, *udot;

    149:   PetscFunctionBegin;
    150:   /*  The next three lines allow us to access the entries of the vectors directly */
    151:   PetscCall(VecGetArrayRead(U, &u));
    152:   PetscCall(VecGetArrayRead(Udot, &udot));
    153:   PetscCall(VecGetArray(F, &f));

    155:   if (actx->mode == 1) {
    156:     f[0] = udot[0] - u[0] + 100 * u[1];
    157:     f[1] = udot[1] - 10 * u[0] - u[1];
    158:   } else if (actx->mode == 2) {
    159:     f[0] = udot[0] - u[0] - 10 * u[1];
    160:     f[1] = udot[1] + 100 * u[0] - u[1];
    161:   }

    163:   PetscCall(VecRestoreArrayRead(U, &u));
    164:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    165:   PetscCall(VecRestoreArray(F, &f));
    166:   PetscFunctionReturn(PETSC_SUCCESS);
    167: }

    169: /*
    170:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
    171: */
    172: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
    173: {
    174:   AppCtx            *actx     = (AppCtx *)ctx;
    175:   PetscInt           rowcol[] = {0, 1};
    176:   PetscScalar        J[2][2];
    177:   const PetscScalar *u, *udot;

    179:   PetscFunctionBegin;
    180:   PetscCall(VecGetArrayRead(U, &u));
    181:   PetscCall(VecGetArrayRead(Udot, &udot));

    183:   if (actx->mode == 1) {
    184:     J[0][0] = a - 1;
    185:     J[0][1] = 100;
    186:     J[1][0] = -10;
    187:     J[1][1] = a - 1;
    188:   } else if (actx->mode == 2) {
    189:     J[0][0] = a - 1;
    190:     J[0][1] = -10;
    191:     J[1][0] = 100;
    192:     J[1][1] = a - 1;
    193:   }
    194:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));

    196:   PetscCall(VecRestoreArrayRead(U, &u));
    197:   PetscCall(VecRestoreArrayRead(Udot, &udot));

    199:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    200:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    201:   if (A != B) {
    202:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    203:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    204:   }
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
    209: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat Ap, void *ctx)
    210: {
    211:   PetscFunctionBeginUser;
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: int main(int argc, char **argv)
    216: {
    217:   TS           ts; /* ODE integrator */
    218:   Vec          U;  /* solution will be stored here */
    219:   Mat          A;  /* Jacobian matrix */
    220:   Mat          Ap; /* Jacobian dfdp */
    221:   PetscMPIInt  size;
    222:   PetscInt     n = 2;
    223:   PetscScalar *u;
    224:   AppCtx       app;
    225:   PetscInt     direction[1];
    226:   PetscBool    terminate[1];
    227:   Mat          sp;
    228:   PetscReal    tend;
    229:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    230:      Initialize program
    231:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    232:   PetscFunctionBeginUser;
    233:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    234:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    235:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
    236:   app.mode       = 1;
    237:   app.lambda1    = 2.75;
    238:   app.lambda2    = 0.36;
    239:   app.print_time = 1. / 256.;
    240:   tend           = 0.125;
    241:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex1fwd options", "");
    242:   {
    243:     PetscCall(PetscOptionsReal("-lambda1", "", "", app.lambda1, &app.lambda1, NULL));
    244:     PetscCall(PetscOptionsReal("-lambda2", "", "", app.lambda2, &app.lambda2, NULL));
    245:     PetscCall(PetscOptionsReal("-tend", "", "", tend, &tend, NULL));
    246:   }
    247:   PetscOptionsEnd();

    249:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    250:     Create necessary matrix and vectors
    251:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    252:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    253:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    254:   PetscCall(MatSetType(A, MATDENSE));
    255:   PetscCall(MatSetFromOptions(A));
    256:   PetscCall(MatSetUp(A));

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

    260:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Ap));
    261:   PetscCall(MatSetSizes(Ap, n, 3, PETSC_DETERMINE, PETSC_DETERMINE));
    262:   PetscCall(MatSetType(Ap, MATDENSE));
    263:   PetscCall(MatSetFromOptions(Ap));
    264:   PetscCall(MatSetUp(Ap));
    265:   PetscCall(MatZeroEntries(Ap));

    267:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, 3, NULL, &sp));
    268:   PetscCall(MatZeroEntries(sp));
    269:   PetscCall(MatShift(sp, 1.0));

    271:   PetscCall(VecGetArray(U, &u));
    272:   u[0] = 0;
    273:   u[1] = 1;
    274:   PetscCall(VecRestoreArray(U, &u));

    276:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    277:      Create timestepping solver context
    278:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    279:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    280:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    281:   PetscCall(TSSetType(ts, TSCN));
    282:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &app));
    283:   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &app));

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

    290:   PetscCall(TSForwardSetSensitivities(ts, 3, sp));
    291:   /*   Set RHS JacobianP */
    292:   PetscCall(TSSetRHSJacobianP(ts, Ap, RHSJacobianP, &app));

    294:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    295:      Set solver options
    296:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    297:   PetscCall(TSSetMaxTime(ts, tend));
    298:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    299:   PetscCall(TSSetTimeStep(ts, 1. / 256.));
    300:   PetscCall(TSMonitorSet(ts, MyMonitor, &app, NULL));
    301:   PetscCall(TSSetFromOptions(ts));

    303:   /* Set directions and terminate flags for the two events */
    304:   direction[0] = 0;
    305:   terminate[0] = PETSC_FALSE;
    306:   PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));
    307:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    308:      Run timestepping solver
    309:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    310:   PetscCall(TSSolve(ts, U));

    312:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    313:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    314:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    315:   PetscCall(MatDestroy(&A));
    316:   PetscCall(VecDestroy(&U));
    317:   PetscCall(TSDestroy(&ts));

    319:   PetscCall(MatDestroy(&Ap));
    320:   PetscCall(MatDestroy(&sp));
    321:   PetscCall(PetscFinalize());
    322:   return 0;
    323: }

    325: /*TEST

    327:    build:
    328:       requires: !complex

    330:    test:
    331:       args: -ts_monitor

    333: TEST*/