Actual source code: drivers.cxx

      1: #include "contexts.cxx"
      2: #include "sparse.cxx"
      3: #include "init.cxx"
      4: #include <adolc/drivers/drivers.h>
      5: #include <adolc/interfaces.h>

      7: /*
      8:    REQUIRES configuration of PETSc with option --download-adolc.

     10:    For documentation on ADOL-C, see
     11:      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
     12: */

     14: /* --------------------------------------------------------------------------------
     15:    Drivers for RHSJacobian and IJacobian
     16:    ----------------------------------------------------------------------------- */

     18: /*
     19:   Compute Jacobian for explicit TS in compressed format and recover from this, using
     20:   precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
     21:   assembled (not recommended for non-toy problems!).

     23:   Input parameters:
     24:   tag   - tape identifier
     25:   u_vec - vector at which to evaluate Jacobian
     26:   ctx   - ADOL-C context, as defined above

     28:   Output parameter:
     29:   A     - Mat object corresponding to Jacobian
     30: */
     31: PetscErrorCode PetscAdolcComputeRHSJacobian(PetscInt tag, Mat A, const PetscScalar *u_vec, void *ctx)
     32: {
     33:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
     34:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
     35:   PetscScalar **J;

     37:   PetscFunctionBegin;
     38:   PetscCall(AdolcMalloc2(m, p, &J));
     39:   if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
     40:   else jacobian(tag, m, n, u_vec, J);
     41:   if (adctx->sparse) {
     42:     PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
     43:   } else {
     44:     for (i = 0; i < m; i++) {
     45:       for (j = 0; j < n; j++) {
     46:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
     47:       }
     48:     }
     49:   }
     50:   PetscCall(AdolcFree2(J));
     51:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     52:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: /*
     57:   Compute Jacobian for explicit TS in compressed format and recover from this, using
     58:   precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
     59:   assembled (not recommended for non-toy problems!).

     61:   Input parameters:
     62:   tag   - tape identifier
     63:   u_vec - vector at which to evaluate Jacobian
     64:   ctx   - ADOL-C context, as defined above

     66:   Output parameter:
     67:   A     - Mat object corresponding to Jacobian
     68: */
     69: PetscErrorCode PetscAdolcComputeRHSJacobianLocal(PetscInt tag, Mat A, const PetscScalar *u_vec, void *ctx)
     70: {
     71:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
     72:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
     73:   PetscScalar **J;

     75:   PetscFunctionBegin;
     76:   PetscCall(AdolcMalloc2(m, p, &J));
     77:   if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
     78:   else jacobian(tag, m, n, u_vec, J);
     79:   if (adctx->sparse) {
     80:     PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
     81:   } else {
     82:     for (i = 0; i < m; i++) {
     83:       for (j = 0; j < n; j++) {
     84:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
     85:       }
     86:     }
     87:   }
     88:   PetscCall(AdolcFree2(J));
     89:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     90:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: /*
     95:   Compute Jacobian for implicit TS in compressed format and recover from this, using
     96:   precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
     97:   assembled (not recommended for non-toy problems!).

     99:   Input parameters:
    100:   tag1   - tape identifier for df/dx part
    101:   tag2   - tape identifier for df/d(xdot) part
    102:   u_vec - vector at which to evaluate Jacobian
    103:   ctx   - ADOL-C context, as defined above

    105:   Output parameter:
    106:   A     - Mat object corresponding to Jacobian
    107: */
    108: PetscErrorCode PetscAdolcComputeIJacobian(PetscInt tag1, PetscInt tag2, Mat A, const PetscScalar *u_vec, PetscReal a, void *ctx)
    109: {
    110:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    111:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
    112:   PetscScalar **J;

    114:   PetscFunctionBegin;
    115:   PetscCall(AdolcMalloc2(m, p, &J));

    117:   /* dF/dx part */
    118:   if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
    119:   else jacobian(tag1, m, n, u_vec, J);
    120:   PetscCall(MatZeroEntries(A));
    121:   if (adctx->sparse) {
    122:     PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
    123:   } else {
    124:     for (i = 0; i < m; i++) {
    125:       for (j = 0; j < n; j++) {
    126:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    127:       }
    128:     }
    129:   }
    130:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    131:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    133:   /* a * dF/d(xdot) part */
    134:   if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
    135:   else jacobian(tag2, m, n, u_vec, J);
    136:   if (adctx->sparse) {
    137:     PetscCall(RecoverJacobian(A, ADD_VALUES, m, p, adctx->Rec, J, &a));
    138:   } else {
    139:     for (i = 0; i < m; i++) {
    140:       for (j = 0; j < n; j++) {
    141:         if (fabs(J[i][j]) > 1.e-16) {
    142:           J[i][j] *= a;
    143:           PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], ADD_VALUES));
    144:         }
    145:       }
    146:     }
    147:   }
    148:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    149:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    150:   PetscCall(AdolcFree2(J));
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: /*
    155:   Compute Jacobian for implicit TS in the special case where it is
    156:   known that the mass matrix is simply the identity. i.e. We have
    157:   a problem of the form
    158:                         du/dt = F(u).

    160:   Input parameters:
    161:   tag   - tape identifier for df/dx part
    162:   u_vec - vector at which to evaluate Jacobian
    163:   ctx   - ADOL-C context, as defined above

    165:   Output parameter:
    166:   A     - Mat object corresponding to Jacobian
    167: */
    168: PetscErrorCode PetscAdolcComputeIJacobianIDMass(PetscInt tag, Mat A, PetscScalar *u_vec, PetscReal a, void *ctx)
    169: {
    170:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    171:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
    172:   PetscScalar **J;

    174:   PetscFunctionBegin;
    175:   PetscCall(AdolcMalloc2(m, p, &J));

    177:   /* dF/dx part */
    178:   if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
    179:   else jacobian(tag, m, n, u_vec, J);
    180:   PetscCall(MatZeroEntries(A));
    181:   if (adctx->sparse) {
    182:     PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
    183:   } else {
    184:     for (i = 0; i < m; i++) {
    185:       for (j = 0; j < n; j++) {
    186:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    187:       }
    188:     }
    189:   }
    190:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    191:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    192:   PetscCall(AdolcFree2(J));

    194:   /* a * dF/d(xdot) part */
    195:   PetscCall(MatShift(A, a));
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: /*
    200:   Compute local portion of Jacobian for implicit TS in compressed format and recover from this, using
    201:   precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
    202:   assembled (not recommended for non-toy problems!).

    204:   Input parameters:
    205:   tag1   - tape identifier for df/dx part
    206:   tag2   - tape identifier for df/d(xdot) part
    207:   u_vec - vector at which to evaluate Jacobian
    208:   ctx   - ADOL-C context, as defined above

    210:   Output parameter:
    211:   A     - Mat object corresponding to Jacobian
    212: */
    213: PetscErrorCode PetscAdolcComputeIJacobianLocal(PetscInt tag1, PetscInt tag2, Mat A, PetscScalar *u_vec, PetscReal a, void *ctx)
    214: {
    215:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    216:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
    217:   PetscScalar **J;

    219:   PetscFunctionBegin;
    220:   PetscCall(AdolcMalloc2(m, p, &J));

    222:   /* dF/dx part */
    223:   if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
    224:   else jacobian(tag1, m, n, u_vec, J);
    225:   if (adctx->sparse) {
    226:     PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
    227:   } else {
    228:     for (i = 0; i < m; i++) {
    229:       for (j = 0; j < n; j++) {
    230:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    231:       }
    232:     }
    233:   }
    234:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    235:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    237:   /* a * dF/d(xdot) part */
    238:   if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
    239:   else jacobian(tag2, m, n, u_vec, J);
    240:   if (adctx->sparse) {
    241:     PetscCall(RecoverJacobianLocal(A, ADD_VALUES, m, p, adctx->Rec, J, &a));
    242:   } else {
    243:     for (i = 0; i < m; i++) {
    244:       for (j = 0; j < n; j++) {
    245:         if (fabs(J[i][j]) > 1.e-16) {
    246:           J[i][j] *= a;
    247:           PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], ADD_VALUES));
    248:         }
    249:       }
    250:     }
    251:   }
    252:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    253:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    254:   PetscCall(AdolcFree2(J));
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: /*
    259:   Compute local portion of Jacobian for implicit TS in the special case where it is
    260:   known that the mass matrix is simply the identity. i.e. We have
    261:   a problem of the form
    262:                         du/dt = F(u).

    264:   Input parameters:
    265:   tag   - tape identifier for df/dx part
    266:   u_vec - vector at which to evaluate Jacobian
    267:   ctx   - ADOL-C context, as defined above

    269:   Output parameter:
    270:   A     - Mat object corresponding to Jacobian
    271: */
    272: PetscErrorCode PetscAdolcComputeIJacobianLocalIDMass(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscReal a, void *ctx)
    273: {
    274:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    275:   PetscInt      i, j, m = adctx->m, n = adctx->n, p = adctx->p;
    276:   PetscScalar **J;

    278:   PetscFunctionBegin;
    279:   PetscCall(AdolcMalloc2(m, p, &J));

    281:   /* dF/dx part */
    282:   if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
    283:   else jacobian(tag, m, n, u_vec, J);
    284:   if (adctx->sparse) {
    285:     PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
    286:   } else {
    287:     for (i = 0; i < m; i++) {
    288:       for (j = 0; j < n; j++) {
    289:         if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    290:       }
    291:     }
    292:   }
    293:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    294:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    295:   PetscCall(AdolcFree2(J));

    297:   /* a * dF/d(xdot) part */
    298:   PetscCall(MatShift(A, a));
    299:   PetscFunctionReturn(PETSC_SUCCESS);
    300: }

    302: /* --------------------------------------------------------------------------------
    303:    Drivers for Jacobian w.r.t. a parameter
    304:    ----------------------------------------------------------------------------- */

    306: /*
    307:   Compute Jacobian w.r.t a parameter for explicit TS.

    309:   Input parameters:
    310:   tag    - tape identifier
    311:   u_vec  - vector at which to evaluate Jacobian
    312:   params - the parameters w.r.t. which we differentiate
    313:   ctx    - ADOL-C context, as defined above

    315:   Output parameter:
    316:   A      - Mat object corresponding to Jacobian
    317: */
    318: PetscErrorCode PetscAdolcComputeRHSJacobianP(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscScalar *params, void *ctx)
    319: {
    320:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    321:   PetscInt      i, j = 0, m = adctx->m, n = adctx->n, p = adctx->num_params;
    322:   PetscScalar **J, *concat, **S;

    324:   PetscFunctionBegin;
    325:   /* Allocate memory and concatenate independent variable values with parameter */
    326:   PetscCall(AdolcMalloc2(m, p, &J));
    327:   PetscCall(PetscMalloc1(n + p, &concat));
    328:   PetscCall(AdolcMalloc2(n + p, p, &S));
    329:   PetscCall(Subidentity(p, n, S));
    330:   for (i = 0; i < n; i++) concat[i] = u_vec[i];
    331:   for (i = 0; i < p; i++) concat[n + i] = params[i];

    333:   /* Propagate the appropriate seed matrix through the forward mode of AD */
    334:   fov_forward(tag, m, n + p, p, concat, S, NULL, J);
    335:   PetscCall(AdolcFree2(S));
    336:   PetscCall(PetscFree(concat));

    338:   /* Set matrix values */
    339:   for (i = 0; i < m; i++) {
    340:     for (j = 0; j < p; j++) {
    341:       if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    342:     }
    343:   }
    344:   PetscCall(AdolcFree2(J));
    345:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    346:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    347:   PetscFunctionReturn(PETSC_SUCCESS);
    348: }

    350: /*
    351:   Compute local portion of Jacobian w.r.t a parameter for explicit TS.

    353:   Input parameters:
    354:   tag    - tape identifier
    355:   u_vec  - vector at which to evaluate Jacobian
    356:   params - the parameters w.r.t. which we differentiate
    357:   ctx    - ADOL-C context, as defined above

    359:   Output parameter:
    360:   A      - Mat object corresponding to Jacobian
    361: */
    362: PetscErrorCode PetscAdolcComputeRHSJacobianPLocal(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscScalar *params, void *ctx)
    363: {
    364:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    365:   PetscInt      i, j = 0, m = adctx->m, n = adctx->n, p = adctx->num_params;
    366:   PetscScalar **J, *concat, **S;

    368:   PetscFunctionBegin;
    369:   /* Allocate memory and concatenate independent variable values with parameter */
    370:   PetscCall(AdolcMalloc2(m, p, &J));
    371:   PetscCall(PetscMalloc1(n + p, &concat));
    372:   PetscCall(AdolcMalloc2(n + p, p, &S));
    373:   PetscCall(Subidentity(p, n, S));
    374:   for (i = 0; i < n; i++) concat[i] = u_vec[i];
    375:   for (i = 0; i < p; i++) concat[n + i] = params[i];

    377:   /* Propagate the appropriate seed matrix through the forward mode of AD */
    378:   fov_forward(tag, m, n + p, p, concat, S, NULL, J);
    379:   PetscCall(AdolcFree2(S));
    380:   PetscCall(PetscFree(concat));

    382:   /* Set matrix values */
    383:   for (i = 0; i < m; i++) {
    384:     for (j = 0; j < p; j++) {
    385:       if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
    386:     }
    387:   }
    388:   PetscCall(AdolcFree2(J));
    389:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    390:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    391:   PetscFunctionReturn(PETSC_SUCCESS);
    392: }

    394: /* --------------------------------------------------------------------------------
    395:    Drivers for Jacobian diagonal
    396:    ----------------------------------------------------------------------------- */

    398: /*
    399:   Compute local portion of Jacobian diagonal for implicit TS in compressed format and recover
    400:   from this, using precomputed seed matrix and recovery vector.

    402:   Input parameters:
    403:   tag1  - tape identifier for df/dx part
    404:   tag2  - tape identifier for df/d(xdot) part
    405:   u_vec - vector at which to evaluate Jacobian
    406:   ctx   - ADOL-C context, as defined above

    408:   Output parameter:
    409:   diag  - Vec object corresponding to Jacobian diagonal
    410: */
    411: PetscErrorCode PetscAdolcComputeIJacobianAndDiagonalLocal(PetscInt tag1, PetscInt tag2, Vec diag, PetscScalar *u_vec, PetscReal a, void *ctx)
    412: {
    413:   AdolcCtx     *adctx = (AdolcCtx *)ctx;
    414:   PetscInt      i, m = adctx->m, n = adctx->n, p = adctx->p;
    415:   PetscScalar **J;

    417:   PetscFunctionBegin;
    418:   PetscCall(AdolcMalloc2(m, p, &J));

    420:   /* dF/dx part */
    421:   if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
    422:   else jacobian(tag1, m, n, u_vec, J);
    423:   if (adctx->sparse) {
    424:     PetscCall(RecoverDiagonalLocal(diag, INSERT_VALUES, m, adctx->rec, J, NULL));
    425:   } else {
    426:     for (i = 0; i < m; i++) {
    427:       if (fabs(J[i][i]) > 1.e-16) PetscCall(VecSetValuesLocal(diag, 1, &i, &J[i][i], INSERT_VALUES));
    428:     }
    429:   }
    430:   PetscCall(VecAssemblyBegin(diag));
    431:   PetscCall(VecAssemblyEnd(diag));

    433:   /* a * dF/d(xdot) part */
    434:   if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
    435:   else jacobian(tag2, m, n, u_vec, J);
    436:   if (adctx->sparse) {
    437:     PetscCall(RecoverDiagonalLocal(diag, ADD_VALUES, m, adctx->rec, J, NULL));
    438:   } else {
    439:     for (i = 0; i < m; i++) {
    440:       if (fabs(J[i][i]) > 1.e-16) {
    441:         J[i][i] *= a;
    442:         PetscCall(VecSetValuesLocal(diag, 1, &i, &J[i][i], ADD_VALUES));
    443:       }
    444:     }
    445:   }
    446:   PetscCall(VecAssemblyBegin(diag));
    447:   PetscCall(VecAssemblyEnd(diag));
    448:   PetscCall(AdolcFree2(J));
    449:   PetscFunctionReturn(PETSC_SUCCESS);
    450: }