Actual source code: adr_ex1.cxx

      1: static char help[] = "Demonstrates automatic Jacobian generation using ADOL-C for a nonlinear reaction problem from chemistry.\n";

      3: /*
      4:    REQUIRES configuration of PETSc with option --download-adolc.

      6:    For documentation on ADOL-C, see
      7:      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
      8: */
      9: /* ------------------------------------------------------------------------
     10:   See ../advection-diffusion-reaction/ex1 for a description of the problem
     11:   ------------------------------------------------------------------------- */
     12: #include <petscts.h>
     13: #include "adolc-utils/drivers.cxx"
     14: #include <adolc/adolc.h>

     16: typedef struct {
     17:   PetscScalar k;
     18:   Vec         initialsolution;
     19:   AdolcCtx   *adctx; /* Automatic differentiation support */
     20: } AppCtx;

     22: PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v)
     23: {
     24:   PetscFunctionBegin;
     25:   PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v)
     30: {
     31:   PetscFunctionBegin;
     32:   PetscCall(PetscNew(ctx));
     33:   PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: /*
     38:   Defines the ODE passed to the ODE solver
     39: */
     40: PetscErrorCode IFunctionPassive(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
     41: {
     42:   PetscScalar       *f;
     43:   const PetscScalar *u, *udot;

     45:   PetscFunctionBegin;
     46:   /*  The next three lines allow us to access the entries of the vectors directly */
     47:   PetscCall(VecGetArrayRead(U, &u));
     48:   PetscCall(VecGetArrayRead(Udot, &udot));
     49:   PetscCall(VecGetArray(F, &f));
     50:   f[0] = udot[0] + ctx->k * u[0] * u[1];
     51:   f[1] = udot[1] + ctx->k * u[0] * u[1];
     52:   f[2] = udot[2] - ctx->k * u[0] * u[1];
     53:   PetscCall(VecRestoreArray(F, &f));
     54:   PetscCall(VecRestoreArrayRead(Udot, &udot));
     55:   PetscCall(VecRestoreArrayRead(U, &u));
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: /*
     60:   'Active' ADOL-C annotated version, marking dependence upon u.
     61: */
     62: PetscErrorCode IFunctionActive1(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
     63: {
     64:   PetscScalar       *f;
     65:   const PetscScalar *u, *udot;

     67:   adouble f_a[3]; /* 'active' double for dependent variables */
     68:   adouble u_a[3]; /* 'active' double for independent variables */

     70:   PetscFunctionBegin;
     71:   /*  The next three lines allow us to access the entries of the vectors directly */
     72:   PetscCall(VecGetArrayRead(U, &u));
     73:   PetscCall(VecGetArrayRead(Udot, &udot));
     74:   PetscCall(VecGetArray(F, &f));

     76:   /* Start of active section */
     77:   trace_on(1);
     78:   u_a[0] <<= u[0];
     79:   u_a[1] <<= u[1];
     80:   u_a[2] <<= u[2]; /* Mark independence */
     81:   f_a[0] = udot[0] + ctx->k * u_a[0] * u_a[1];
     82:   f_a[1] = udot[1] + ctx->k * u_a[0] * u_a[1];
     83:   f_a[2] = udot[2] - ctx->k * u_a[0] * u_a[1];
     84:   f_a[0] >>= f[0];
     85:   f_a[1] >>= f[1];
     86:   f_a[2] >>= f[2]; /* Mark dependence */
     87:   trace_off();
     88:   /* End of active section */

     90:   PetscCall(VecRestoreArray(F, &f));
     91:   PetscCall(VecRestoreArrayRead(Udot, &udot));
     92:   PetscCall(VecRestoreArrayRead(U, &u));
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: /*
     97:   'Active' ADOL-C annotated version, marking dependence upon udot.
     98: */
     99: PetscErrorCode IFunctionActive2(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
    100: {
    101:   PetscScalar       *f;
    102:   const PetscScalar *u, *udot;

    104:   adouble f_a[3];    /* 'active' double for dependent variables */
    105:   adouble udot_a[3]; /* 'active' double for independent variables */

    107:   PetscFunctionBegin;
    108:   /*  The next three lines allow us to access the entries of the vectors directly */
    109:   PetscCall(VecGetArrayRead(U, &u));
    110:   PetscCall(VecGetArrayRead(Udot, &udot));
    111:   PetscCall(VecGetArray(F, &f));

    113:   /* Start of active section */
    114:   trace_on(2);
    115:   udot_a[0] <<= udot[0];
    116:   udot_a[1] <<= udot[1];
    117:   udot_a[2] <<= udot[2]; /* Mark independence */
    118:   f_a[0] = udot_a[0] + ctx->k * u[0] * u[1];
    119:   f_a[1] = udot_a[1] + ctx->k * u[0] * u[1];
    120:   f_a[2] = udot_a[2] - ctx->k * u[0] * u[1];
    121:   f_a[0] >>= f[0];
    122:   f_a[1] >>= f[1];
    123:   f_a[2] >>= f[2]; /* Mark dependence */
    124:   trace_off();
    125:   /* End of active section */

    127:   PetscCall(VecRestoreArray(F, &f));
    128:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    129:   PetscCall(VecRestoreArrayRead(U, &u));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*
    134:  Defines the Jacobian of the ODE passed to the ODE solver, using the PETSc-ADOL-C driver for
    135:  implicit TS.
    136: */
    137: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
    138: {
    139:   AppCtx            *appctx = (AppCtx *)ctx;
    140:   const PetscScalar *u;

    142:   PetscFunctionBegin;
    143:   PetscCall(VecGetArrayRead(U, &u));
    144:   PetscCall(PetscAdolcComputeIJacobian(1, 2, A, u, a, appctx->adctx));
    145:   PetscCall(VecRestoreArrayRead(U, &u));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*
    150:      Defines the exact (analytic) solution to the ODE
    151: */
    152: static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
    153: {
    154:   const PetscScalar *uinit;
    155:   PetscScalar       *u, d0, q;

    157:   PetscFunctionBegin;
    158:   PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit));
    159:   PetscCall(VecGetArray(U, &u));
    160:   d0 = uinit[0] - uinit[1];
    161:   if (d0 == 0.0) q = ctx->k * t;
    162:   else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0;
    163:   u[0] = uinit[0] / (1.0 + uinit[1] * q);
    164:   u[1] = u[0] - d0;
    165:   u[2] = uinit[1] + uinit[2] - u[1];
    166:   PetscCall(VecRestoreArray(U, &u));
    167:   PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit));
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: int main(int argc, char **argv)
    172: {
    173:   TS                ts;         /* ODE integrator */
    174:   Vec               U, Udot, R; /* solution, derivative, residual */
    175:   Mat               A;          /* Jacobian matrix */
    176:   PetscMPIInt       size;
    177:   PetscInt          n = 3;
    178:   AppCtx            ctx;
    179:   AdolcCtx         *adctx;
    180:   PetscScalar      *u;
    181:   const char *const names[] = {"U1", "U2", "U3", NULL};

    183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    184:      Initialize program
    185:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    186:   PetscFunctionBeginUser;
    187:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    188:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    189:   PetscCheck(size <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Only for sequential runs");
    190:   PetscCall(PetscNew(&adctx));
    191:   adctx->m  = n;
    192:   adctx->n  = n;
    193:   adctx->p  = n;
    194:   ctx.adctx = adctx;

    196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    197:     Create necessary matrix and vectors
    198:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    199:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    200:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    201:   PetscCall(MatSetFromOptions(A));
    202:   PetscCall(MatSetUp(A));

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

    206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    207:     Set runtime options
    208:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    209:   ctx.k = .9;
    210:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL));
    211:   PetscCall(VecDuplicate(U, &ctx.initialsolution));
    212:   PetscCall(VecGetArray(ctx.initialsolution, &u));
    213:   u[0] = 1;
    214:   u[1] = .7;
    215:   u[2] = 0;
    216:   PetscCall(VecRestoreArray(ctx.initialsolution, &u));
    217:   PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL));

    219:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    220:      Create timestepping solver context
    221:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    222:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    223:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    224:   PetscCall(TSSetType(ts, TSROSW));
    225:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunctionPassive, &ctx));

    227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    228:      Set initial conditions
    229:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    230:   PetscCall(Solution(ts, 0, U, &ctx));
    231:   PetscCall(TSSetSolution(ts, U));

    233:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    234:      Trace just once for each tape
    235:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    236:   PetscCall(VecDuplicate(U, &Udot));
    237:   PetscCall(VecDuplicate(U, &R));
    238:   PetscCall(IFunctionActive1(ts, 0., U, Udot, R, &ctx));
    239:   PetscCall(IFunctionActive2(ts, 0., U, Udot, R, &ctx));
    240:   PetscCall(VecDestroy(&R));
    241:   PetscCall(VecDestroy(&Udot));

    243:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    244:      Set Jacobian
    245:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    246:   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));
    247:   PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx));

    249:   {
    250:     DM    dm;
    251:     void *ptr;
    252:     PetscCall(TSGetDM(ts, &dm));
    253:     PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr));
    254:     PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr));
    255:     PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
    256:     PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
    257:   }

    259:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    260:      Set solver options
    261:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    262:   PetscCall(TSSetTimeStep(ts, .001));
    263:   PetscCall(TSSetMaxSteps(ts, 1000));
    264:   PetscCall(TSSetMaxTime(ts, 20.0));
    265:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    266:   PetscCall(TSSetFromOptions(ts));
    267:   PetscCall(TSMonitorLGSetVariableNames(ts, names));

    269:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    270:      Solve nonlinear system
    271:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    272:   PetscCall(TSSolve(ts, U));

    274:   PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));

    276:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    277:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    278:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    279:   PetscCall(VecDestroy(&ctx.initialsolution));
    280:   PetscCall(MatDestroy(&A));
    281:   PetscCall(VecDestroy(&U));
    282:   PetscCall(TSDestroy(&ts));
    283:   PetscCall(PetscFree(adctx));
    284:   PetscCall(PetscFinalize());
    285:   return 0;
    286: }

    288: /*TEST

    290:   build:
    291:     requires: double !complex adolc

    293:   test:
    294:     suffix: 1
    295:     args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor
    296:     output_file: output/adr_ex1_1.out

    298:   test:
    299:     suffix: 2
    300:     args: -ts_max_steps 1 -snes_test_jacobian
    301:     output_file: output/adr_ex1_2.out

    303: TEST*/