Actual source code: ex16adj_tl.cxx

      1: static char help[] = "Demonstrates tapeless automatic Jacobian generation using ADOL-C for an adjoint sensitivity analysis of the van der Pol equation.\n\
      2: Input parameters include:\n\
      3:       -mu : stiffness parameter\n\n";

      5: /*
      6:    REQUIRES configuration of PETSc with option --download-adolc.

      8:    For documentation on ADOL-C, see
      9:      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
     10: */
     11: /* ------------------------------------------------------------------------
     12:    See ex16adj for a description of the problem being solved.
     13:   ------------------------------------------------------------------------- */

     15: #include <petscts.h>
     16: #include <petscmat.h>

     18: #define ADOLC_TAPELESS
     19: #define NUMBER_DIRECTIONS 3
     20: #include "adolc-utils/drivers.cxx"
     21: #include <adolc/adtl.h>
     22: using namespace adtl;

     24: typedef struct _n_User *User;
     25: struct _n_User {
     26:   PetscReal mu;
     27:   PetscReal next_output;
     28:   PetscReal tprev;

     30:   /* Automatic differentiation support */
     31:   AdolcCtx *adctx;
     32:   Vec       F;
     33: };

     35: /*
     36:   Residual evaluation templated, so as to allow for PetscScalar or adouble
     37:   arguments.
     38: */
     39: template <class T>
     40: PetscErrorCode EvaluateResidual(const T *x, T mu, T *f)
     41: {
     42:   PetscFunctionBegin;
     43:   f[0] = x[1];
     44:   f[1] = mu * (1. - x[0] * x[0]) * x[1] - x[0];
     45:   PetscFunctionReturn(PETSC_SUCCESS);
     46: }

     48: /*
     49:   'Passive' RHS function, used in residual evaluations during the time integration.
     50: */
     51: static PetscErrorCode RHSFunctionPassive(TS ts, PetscReal t, Vec X, Vec F, PetscCtx ctx)
     52: {
     53:   User               user = (User)ctx;
     54:   PetscScalar       *f;
     55:   const PetscScalar *x;

     57:   PetscFunctionBeginUser;
     58:   PetscCall(VecGetArrayRead(X, &x));
     59:   PetscCall(VecGetArray(F, &f));
     60:   PetscCall(EvaluateResidual(x, user->mu, f));
     61:   PetscCall(VecRestoreArrayRead(X, &x));
     62:   PetscCall(VecRestoreArray(F, &f));
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: /*
     67:   Compute the Jacobian w.r.t. x using tapeless mode of ADOL-C.
     68: */
     69: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat A, Mat B, PetscCtx ctx)
     70: {
     71:   User               user = (User)ctx;
     72:   PetscScalar      **J;
     73:   const PetscScalar *x;
     74:   adouble            f_a[2];       /* 'active' double for dependent variables */
     75:   adouble            x_a[2], mu_a; /* 'active' doubles for independent variables */
     76:   PetscInt           i, j;

     78:   PetscFunctionBeginUser;
     79:   /* Set values for independent variables and parameters */
     80:   PetscCall(VecGetArrayRead(X, &x));
     81:   x_a[0].setValue(x[0]);
     82:   x_a[1].setValue(x[1]);
     83:   mu_a.setValue(user->mu);
     84:   PetscCall(VecRestoreArrayRead(X, &x));

     86:   /* Set seed matrix as 3x3 identity matrix */
     87:   x_a[0].setADValue(0, 1.);
     88:   x_a[0].setADValue(1, 0.);
     89:   x_a[0].setADValue(2, 0.);
     90:   x_a[1].setADValue(0, 0.);
     91:   x_a[1].setADValue(1, 1.);
     92:   x_a[1].setADValue(2, 0.);
     93:   mu_a.setADValue(0, 0.);
     94:   mu_a.setADValue(1, 0.);
     95:   mu_a.setADValue(2, 1.);

     97:   /* Evaluate residual (on active variables) */
     98:   PetscCall(EvaluateResidual(x_a, mu_a, f_a));

    100:   /* Extract derivatives */
    101:   PetscCall(PetscMalloc1(user->adctx->n, &J));
    102:   J[0] = (PetscScalar *)f_a[0].getADValue();
    103:   J[1] = (PetscScalar *)f_a[1].getADValue();

    105:   /* Set matrix values */
    106:   for (i = 0; i < user->adctx->m; i++) {
    107:     for (j = 0; j < user->adctx->n; j++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    108:   }
    109:   PetscCall(PetscFree(J));
    110:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    111:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    112:   if (A != B) {
    113:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    114:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    115:   }
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: /*
    120:   Compute the Jacobian w.r.t. mu using tapeless mode of ADOL-C.
    121: */
    122: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx)
    123: {
    124:   User          user = (User)ctx;
    125:   PetscScalar **J;
    126:   PetscScalar  *x;
    127:   adouble       f_a[2];       /* 'active' double for dependent variables */
    128:   adouble       x_a[2], mu_a; /* 'active' doubles for independent variables */
    129:   PetscInt      i, j = 0;

    131:   PetscFunctionBeginUser;
    132:   /* Set values for independent variables and parameters */
    133:   PetscCall(VecGetArray(X, &x));
    134:   x_a[0].setValue(x[0]);
    135:   x_a[1].setValue(x[1]);
    136:   mu_a.setValue(user->mu);
    137:   PetscCall(VecRestoreArray(X, &x));

    139:   /* Set seed matrix as 3x3 identity matrix */
    140:   x_a[0].setADValue(0, 1.);
    141:   x_a[0].setADValue(1, 0.);
    142:   x_a[0].setADValue(2, 0.);
    143:   x_a[1].setADValue(0, 0.);
    144:   x_a[1].setADValue(1, 1.);
    145:   x_a[1].setADValue(2, 0.);
    146:   mu_a.setADValue(0, 0.);
    147:   mu_a.setADValue(1, 0.);
    148:   mu_a.setADValue(2, 1.);

    150:   /* Evaluate residual (on active variables) */
    151:   PetscCall(EvaluateResidual(x_a, mu_a, f_a));

    153:   /* Extract derivatives */
    154:   PetscCall(PetscMalloc1(2, &J));
    155:   J[0] = (PetscScalar *)f_a[0].getADValue();
    156:   J[1] = (PetscScalar *)f_a[1].getADValue();

    158:   /* Set matrix values */
    159:   for (i = 0; i < user->adctx->m; i++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][user->adctx->n], INSERT_VALUES));
    160:   PetscCall(PetscFree(J));
    161:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    162:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }

    166: /*
    167:   Monitor timesteps and use interpolation to output at integer multiples of 0.1
    168: */
    169: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, PetscCtx ctx)
    170: {
    171:   const PetscScalar *x;
    172:   PetscReal          tfinal, dt, tprev;
    173:   User               user = (User)ctx;

    175:   PetscFunctionBeginUser;
    176:   PetscCall(TSGetTimeStep(ts, &dt));
    177:   PetscCall(TSGetMaxTime(ts, &tfinal));
    178:   PetscCall(TSGetPrevTime(ts, &tprev));
    179:   PetscCall(VecGetArrayRead(X, &x));
    180:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
    181:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "t %.6f (tprev = %.6f) \n", (double)t, (double)tprev));
    182:   PetscCall(VecRestoreArrayRead(X, &x));
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: int main(int argc, char **argv)
    187: {
    188:   TS             ts;   /* nonlinear solver */
    189:   Vec            x;    /* solution, residual vectors */
    190:   Mat            A;    /* Jacobian matrix */
    191:   Mat            Jacp; /* JacobianP matrix */
    192:   PetscInt       steps;
    193:   PetscReal      ftime   = 0.5;
    194:   PetscBool      monitor = PETSC_FALSE;
    195:   PetscScalar   *x_ptr;
    196:   PetscMPIInt    size;
    197:   struct _n_User user;
    198:   AdolcCtx      *adctx;
    199:   Vec            lambda[2], mu[2];

    201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    202:      Initialize program
    203:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    204:   PetscFunctionBeginUser;
    205:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    206:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    207:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

    209:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    210:     Set runtime options and create AdolcCtx
    211:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    212:   PetscCall(PetscNew(&adctx));
    213:   user.mu          = 1;
    214:   user.next_output = 0.0;
    215:   adctx->m         = 2;
    216:   adctx->n         = 2;
    217:   adctx->p         = 2;
    218:   user.adctx       = adctx;
    219:   adtl::setNumDir(adctx->n + 1); /* #indep. variables, plus parameters */

    221:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
    222:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));

    224:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    225:     Create necessary matrix and vectors, solve same ODE on every process
    226:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    227:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    228:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
    229:   PetscCall(MatSetFromOptions(A));
    230:   PetscCall(MatSetUp(A));
    231:   PetscCall(MatCreateVecs(A, &x, NULL));

    233:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
    234:   PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
    235:   PetscCall(MatSetFromOptions(Jacp));
    236:   PetscCall(MatSetUp(Jacp));

    238:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    239:      Create timestepping solver context
    240:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    241:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    242:   PetscCall(TSSetType(ts, TSRK));
    243:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionPassive, &user));

    245:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    246:      Set initial conditions
    247:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    248:   PetscCall(VecGetArray(x, &x_ptr));
    249:   x_ptr[0] = 2;
    250:   x_ptr[1] = 0.66666654321;
    251:   PetscCall(VecRestoreArray(x, &x_ptr));

    253:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    254:      Set RHS Jacobian for the adjoint integration
    255:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    256:   PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, &user));
    257:   PetscCall(TSSetMaxTime(ts, ftime));
    258:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    259:   if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
    260:   PetscCall(TSSetTimeStep(ts, .001));

    262:   /*
    263:     Have the TS save its trajectory so that TSAdjointSolve() may be used
    264:   */
    265:   PetscCall(TSSetSaveTrajectory(ts));

    267:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    268:      Set runtime options
    269:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    270:   PetscCall(TSSetFromOptions(ts));

    272:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    273:      Solve nonlinear system
    274:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    275:   PetscCall(TSSolve(ts, x));
    276:   PetscCall(TSGetSolveTime(ts, &ftime));
    277:   PetscCall(TSGetStepNumber(ts, &steps));
    278:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, steps, (double)ftime));
    279:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

    281:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    282:      Start the Adjoint model
    283:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    284:   PetscCall(MatCreateVecs(A, &lambda[0], NULL));
    285:   PetscCall(MatCreateVecs(A, &lambda[1], NULL));
    286:   /*   Reset initial conditions for the adjoint integration */
    287:   PetscCall(VecGetArray(lambda[0], &x_ptr));
    288:   x_ptr[0] = 1.0;
    289:   x_ptr[1] = 0.0;
    290:   PetscCall(VecRestoreArray(lambda[0], &x_ptr));
    291:   PetscCall(VecGetArray(lambda[1], &x_ptr));
    292:   x_ptr[0] = 0.0;
    293:   x_ptr[1] = 1.0;
    294:   PetscCall(VecRestoreArray(lambda[1], &x_ptr));

    296:   PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
    297:   PetscCall(MatCreateVecs(Jacp, &mu[1], NULL));
    298:   PetscCall(VecGetArray(mu[0], &x_ptr));
    299:   x_ptr[0] = 0.0;
    300:   PetscCall(VecRestoreArray(mu[0], &x_ptr));
    301:   PetscCall(VecGetArray(mu[1], &x_ptr));
    302:   x_ptr[0] = 0.0;
    303:   PetscCall(VecRestoreArray(mu[1], &x_ptr));
    304:   PetscCall(TSSetCostGradients(ts, 2, lambda, mu));

    306:   /*   Set RHS JacobianP */
    307:   PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &user));

    309:   PetscCall(TSAdjointSolve(ts));

    311:   PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
    312:   PetscCall(VecView(lambda[1], PETSC_VIEWER_STDOUT_WORLD));
    313:   PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
    314:   PetscCall(VecView(mu[1], PETSC_VIEWER_STDOUT_WORLD));

    316:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    317:      Free work space.  All PETSc objects should be destroyed when they
    318:      are no longer needed.
    319:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    320:   PetscCall(MatDestroy(&A));
    321:   PetscCall(MatDestroy(&Jacp));
    322:   PetscCall(VecDestroy(&x));
    323:   PetscCall(VecDestroy(&lambda[0]));
    324:   PetscCall(VecDestroy(&lambda[1]));
    325:   PetscCall(VecDestroy(&mu[0]));
    326:   PetscCall(VecDestroy(&mu[1]));
    327:   PetscCall(TSDestroy(&ts));
    328:   PetscCall(PetscFree(adctx));
    329:   PetscCall(PetscFinalize());
    330:   return 0;
    331: }

    333: /*TEST

    335:   build:
    336:     requires: double !complex adolc

    338:   test:
    339:     suffix: 1
    340:     args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor
    341:     output_file: output/ex16adj_tl_1.out

    343:   test:
    344:     suffix: 2
    345:     args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -mu 5
    346:     output_file: output/ex16adj_tl_2.out

    348:   test:
    349:     suffix: 3
    350:     args: -ts_max_steps 10 -monitor
    351:     output_file: output/ex16adj_tl_3.out

    353:   test:
    354:     suffix: 4
    355:     args: -ts_max_steps 10 -monitor -mu 5
    356:     output_file: output/ex16adj_tl_4.out

    358: TEST*/