Actual source code: ex1fd.c

      1: static char help[] = "An example of hybrid system using TS event.\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:   Reference:
     15:   I. A. Hiskens, M.A. Pai, Trajectory Sensitivity Analysis of Hybrid Systems, IEEE Transactions on Circuits and Systems, Vol 47, No 2, February 2000
     16: */

     18: #include <petscts.h>

     20: typedef struct {
     21:   PetscReal lambda1;
     22:   PetscReal lambda2;
     23:   PetscInt  mode; /* mode flag*/
     24: } AppCtx;

     26: PetscErrorCode FWDRun(TS, Vec, void *);

     28: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
     29: {
     30:   AppCtx            *actx = (AppCtx *)ctx;
     31:   const PetscScalar *u;

     33:   PetscFunctionBegin;
     34:   PetscCall(VecGetArrayRead(U, &u));
     35:   if (actx->mode == 1) {
     36:     fvalue[0] = PetscRealPart(u[1] - actx->lambda1 * u[0]);
     37:   } else if (actx->mode == 2) {
     38:     fvalue[0] = PetscRealPart(u[1] - actx->lambda2 * u[0]);
     39:   }
     40:   PetscCall(VecRestoreArrayRead(U, &u));
     41:   PetscFunctionReturn(PETSC_SUCCESS);
     42: }

     44: PetscErrorCode ShiftGradients(TS ts, Vec U, AppCtx *actx)
     45: {
     46:   Vec               *lambda, *mu;
     47:   PetscScalar       *x, *y;
     48:   const PetscScalar *u;
     49:   PetscScalar        tmp[2], A1[2][2], A2[2], denorm1, denorm2;
     50:   PetscInt           numcost;

     52:   PetscFunctionBegin;
     53:   PetscCall(TSGetCostGradients(ts, &numcost, &lambda, &mu));
     54:   PetscCall(VecGetArrayRead(U, &u));

     56:   if (actx->mode == 2) {
     57:     denorm1  = -actx->lambda1 * (u[0] - 100. * u[1]) + 1. * (10. * u[0] + u[1]);
     58:     denorm2  = -actx->lambda1 * (u[0] + 10. * u[1]) + 1. * (-100. * u[0] + u[1]);
     59:     A1[0][0] = 110. * u[1] * (-actx->lambda1) / denorm1 + 1.;
     60:     A1[0][1] = -110. * u[0] * (-actx->lambda1) / denorm1;
     61:     A1[1][0] = 110. * u[1] * 1. / denorm1;
     62:     A1[1][1] = -110. * u[0] * 1. / denorm1 + 1.;

     64:     A2[0] = 110. * u[1] * (-u[0]) / denorm2;
     65:     A2[1] = -110. * u[0] * (-u[0]) / denorm2;
     66:   } else {
     67:     denorm2  = -actx->lambda2 * (u[0] + 10. * u[1]) + 1. * (-100. * u[0] + u[1]);
     68:     A1[0][0] = 110. * u[1] * (-actx->lambda1) / denorm2 + 1;
     69:     A1[0][1] = -110. * u[0] * (-actx->lambda1) / denorm2;
     70:     A1[1][0] = 110. * u[1] * 1. / denorm2;
     71:     A1[1][1] = -110. * u[0] * 1. / denorm2 + 1.;

     73:     A2[0] = 0;
     74:     A2[1] = 0;
     75:   }

     77:   PetscCall(VecRestoreArrayRead(U, &u));

     79:   PetscCall(VecGetArray(lambda[0], &x));
     80:   PetscCall(VecGetArray(mu[0], &y));
     81:   tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
     82:   tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
     83:   y[0]   = y[0] + A2[0] * x[0] + A2[1] * x[1];
     84:   x[0]   = tmp[0];
     85:   x[1]   = tmp[1];
     86:   PetscCall(VecRestoreArray(mu[0], &y));
     87:   PetscCall(VecRestoreArray(lambda[0], &x));

     89:   PetscCall(VecGetArray(lambda[1], &x));
     90:   PetscCall(VecGetArray(mu[1], &y));
     91:   tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
     92:   tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
     93:   y[0]   = y[0] + A2[0] * x[0] + A2[1] * x[1];
     94:   x[0]   = tmp[0];
     95:   x[1]   = tmp[1];
     96:   PetscCall(VecRestoreArray(mu[1], &y));
     97:   PetscCall(VecRestoreArray(lambda[1], &x));
     98:   PetscFunctionReturn(PETSC_SUCCESS);
     99: }

    101: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
    102: {
    103:   AppCtx *actx = (AppCtx *)ctx;

    105:   PetscFunctionBegin;
    106:   if (!forwardsolve) PetscCall(ShiftGradients(ts, U, actx));
    107:   if (actx->mode == 1) {
    108:     actx->mode = 2;
    109:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Change from mode 1 to 2 at t = %f \n", (double)t));
    110:   } else if (actx->mode == 2) {
    111:     actx->mode = 1;
    112:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Change from mode 2 to 1 at t = %f \n", (double)t));
    113:   }
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: /*
    118:      Defines the ODE passed to the ODE solver
    119: */
    120: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
    121: {
    122:   AppCtx            *actx = (AppCtx *)ctx;
    123:   PetscScalar       *f;
    124:   const PetscScalar *u, *udot;

    126:   PetscFunctionBegin;
    127:   /*  The next three lines allow us to access the entries of the vectors directly */
    128:   PetscCall(VecGetArrayRead(U, &u));
    129:   PetscCall(VecGetArrayRead(Udot, &udot));
    130:   PetscCall(VecGetArray(F, &f));

    132:   if (actx->mode == 1) {
    133:     f[0] = udot[0] - u[0] + 100 * u[1];
    134:     f[1] = udot[1] - 10 * u[0] - u[1];
    135:   } else if (actx->mode == 2) {
    136:     f[0] = udot[0] - u[0] - 10 * u[1];
    137:     f[1] = udot[1] + 100 * u[0] - u[1];
    138:   }

    140:   PetscCall(VecRestoreArrayRead(U, &u));
    141:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    142:   PetscCall(VecRestoreArray(F, &f));
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

    146: /*
    147:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
    148: */
    149: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
    150: {
    151:   AppCtx            *actx     = (AppCtx *)ctx;
    152:   PetscInt           rowcol[] = {0, 1};
    153:   PetscScalar        J[2][2];
    154:   const PetscScalar *u, *udot;

    156:   PetscFunctionBegin;
    157:   PetscCall(VecGetArrayRead(U, &u));
    158:   PetscCall(VecGetArrayRead(Udot, &udot));

    160:   if (actx->mode == 1) {
    161:     J[0][0] = a - 1;
    162:     J[0][1] = 100;
    163:     J[1][0] = -10;
    164:     J[1][1] = a - 1;
    165:   } else if (actx->mode == 2) {
    166:     J[0][0] = a - 1;
    167:     J[0][1] = -10;
    168:     J[1][0] = 100;
    169:     J[1][1] = a - 1;
    170:   }
    171:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));

    173:   PetscCall(VecRestoreArrayRead(U, &u));
    174:   PetscCall(VecRestoreArrayRead(Udot, &udot));

    176:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    177:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    178:   if (A != B) {
    179:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    180:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    181:   }
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: int main(int argc, char **argv)
    186: {
    187:   TS           ts; /* ODE integrator */
    188:   Vec          U;  /* solution will be stored here */
    189:   Mat          A;  /* Jacobian matrix */
    190:   Mat          Ap; /* dfdp */
    191:   PetscMPIInt  size;
    192:   PetscInt     n = 2;
    193:   PetscScalar *u;
    194:   AppCtx       app;
    195:   PetscInt     direction[1];
    196:   PetscBool    terminate[1];
    197:   PetscReal    delta;
    198:   PetscScalar  tmp[2], sensi[2];

    200:   delta = 1e-8;

    202:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    203:      Initialize program
    204:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    205:   PetscFunctionBeginUser;
    206:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    207:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    208:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
    209:   app.mode    = 1;
    210:   app.lambda1 = 2.75;
    211:   app.lambda2 = 0.36;
    212:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex1 options", "");
    213:   {
    214:     PetscCall(PetscOptionsReal("-lambda1", "", "", app.lambda1, &app.lambda1, NULL));
    215:     PetscCall(PetscOptionsReal("-lambda2", "", "", app.lambda2, &app.lambda2, NULL));
    216:   }
    217:   PetscOptionsEnd();

    219:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    220:     Create necessary matrix and vectors
    221:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    222:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    223:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    224:   PetscCall(MatSetType(A, MATDENSE));
    225:   PetscCall(MatSetFromOptions(A));
    226:   PetscCall(MatSetUp(A));

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

    230:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Ap));
    231:   PetscCall(MatSetSizes(Ap, n, 1, PETSC_DETERMINE, PETSC_DETERMINE));
    232:   PetscCall(MatSetType(Ap, MATDENSE));
    233:   PetscCall(MatSetFromOptions(Ap));
    234:   PetscCall(MatSetUp(Ap));
    235:   PetscCall(MatZeroEntries(Ap)); /* initialize to zeros */

    237:   PetscCall(VecGetArray(U, &u));
    238:   u[0] = 0;
    239:   u[1] = 1;
    240:   PetscCall(VecRestoreArray(U, &u));
    241:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    242:      Create timestepping solver context
    243:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    244:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    245:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    246:   PetscCall(TSSetType(ts, TSCN));
    247:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &app));
    248:   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &app));

    250:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    251:      Set initial conditions
    252:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    253:   PetscCall(TSSetSolution(ts, U));

    255:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    256:      Set solver options
    257:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    258:   PetscCall(TSSetMaxTime(ts, 0.125));
    259:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    260:   PetscCall(TSSetTimeStep(ts, 1. / 256.));
    261:   PetscCall(TSSetFromOptions(ts));

    263:   /* Set directions and terminate flags for the two events */
    264:   direction[0] = 0;
    265:   terminate[0] = PETSC_FALSE;
    266:   PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));

    268:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    269:      Run timestepping solver
    270:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    271:   PetscCall(TSSolve(ts, U));

    273:   PetscCall(VecGetArray(U, &u));
    274:   tmp[0] = u[0];
    275:   tmp[1] = u[1];

    277:   u[0] = 0 + delta;
    278:   u[1] = 1;
    279:   PetscCall(VecRestoreArray(U, &u));

    281:   PetscCall(FWDRun(ts, U, (void *)&app));

    283:   PetscCall(VecGetArray(U, &u));
    284:   sensi[0] = (u[0] - tmp[0]) / delta;
    285:   sensi[1] = (u[1] - tmp[1]) / delta;
    286:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "d x1(tf) /d x1(t0) = %f d x2(tf) / d x1(t0) = %f \n", (double)sensi[0], (double)sensi[1]));
    287:   u[0] = 0;
    288:   u[1] = 1 + delta;
    289:   PetscCall(VecRestoreArray(U, &u));

    291:   PetscCall(FWDRun(ts, U, (void *)&app));

    293:   PetscCall(VecGetArray(U, &u));
    294:   sensi[0] = (u[0] - tmp[0]) / delta;
    295:   sensi[1] = (u[1] - tmp[1]) / delta;
    296:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "d x1(tf) /d x2(t0) = %f d x2(tf) / d x2(t0) = %f \n", (double)sensi[0], (double)sensi[1]));
    297:   u[0]        = 0;
    298:   u[1]        = 1;
    299:   app.lambda1 = app.lambda1 + delta;
    300:   PetscCall(VecRestoreArray(U, &u));

    302:   PetscCall(FWDRun(ts, U, (void *)&app));

    304:   PetscCall(VecGetArray(U, &u));
    305:   sensi[0] = (u[0] - tmp[0]) / delta;
    306:   sensi[1] = (u[1] - tmp[1]) / delta;
    307:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Final gradients: d x1(tf) /d p = %f d x2(tf) / d p = %f \n", (double)sensi[0], (double)sensi[1]));
    308:   PetscCall(VecRestoreArray(U, &u));

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

    317:   PetscCall(MatDestroy(&Ap));
    318:   PetscCall(PetscFinalize());
    319:   return 0;
    320: }

    322: PetscErrorCode FWDRun(TS ts, Vec U0, PetscCtx ctx0)
    323: {
    324:   Vec     U; /* solution will be stored here */
    325:   AppCtx *ctx = (AppCtx *)ctx0;

    327:   PetscFunctionBeginUser;
    328:   PetscCall(TSGetSolution(ts, &U));
    329:   PetscCall(VecCopy(U0, U));

    331:   ctx->mode = 1;
    332:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    333:      Run timestepping solver
    334:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    335:   PetscCall(TSSetTime(ts, 0.0));

    337:   PetscCall(TSSolve(ts, U));
    338:   PetscFunctionReturn(PETSC_SUCCESS);
    339: }

    341: /*TEST

    343:    build:
    344:       requires: !defined(PETSC_USE_CXXCOMPLEX)

    346:    test:
    347:       args: -ts_event_tol 1e-9
    348:       timeoutfactor: 18
    349:       requires: !single

    351: TEST*/