Actual source code: dmadapt.c

      1: #include <petscdmadaptor.h>
      2: #include <petscdmplex.h>
      3: #include <petscdmforest.h>
      4: #include <petscds.h>
      5: #include <petscblaslapack.h>

      7: #include <petsc/private/dmadaptorimpl.h>
      8: #include <petsc/private/dmpleximpl.h>
      9: #include <petsc/private/petscfeimpl.h>

     11: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor, PetscInt, PetscInt, const PetscScalar *, const PetscScalar *, const PetscFVCellGeom *, PetscReal *, void *);

     13: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
     14: {

     18:   DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);
     19:   return(0);
     20: }

     22: /*@
     23:   DMAdaptorCreate - Create a DMAdaptor object. Its purpose is to construct a adaptation DMLabel or metric Vec that can be used to modify the DM.

     25:   Collective

     27:   Input Parameter:
     28: . comm - The communicator for the DMAdaptor object

     30:   Output Parameter:
     31: . adaptor   - The DMAdaptor object

     33:   Level: beginner

     35: .seealso: DMAdaptorDestroy(), DMAdaptorAdapt()
     36: @*/
     37: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
     38: {
     39:   VecTaggerBox   refineBox, coarsenBox;

     44:   PetscSysInitializePackage();
     45:   PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);

     47:   (*adaptor)->monitor = PETSC_FALSE;
     48:   (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
     49:   (*adaptor)->numSeq  = 1;
     50:   (*adaptor)->Nadapt  = -1;
     51:   (*adaptor)->refinementFactor = 2.0;
     52:   (*adaptor)->h_min = 1.;
     53:   (*adaptor)->h_max = 10000.;
     54:   (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
     55:   refineBox.min = refineBox.max = PETSC_MAX_REAL;
     56:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);
     57:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");
     58:   VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);
     59:   VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);
     60:   coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
     61:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);
     62:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");
     63:   VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);
     64:   VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);
     65:   return(0);
     66: }

     68: /*@
     69:   DMAdaptorDestroy - Destroys a DMAdaptor object

     71:   Collective on DMAdaptor

     73:   Input Parameter:
     74: . adaptor - The DMAdaptor object

     76:   Level: beginner

     78: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
     79: @*/
     80: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
     81: {

     85:   if (!*adaptor) return(0);
     87:   if (--((PetscObject)(*adaptor))->refct > 0) {
     88:     *adaptor = NULL;
     89:     return(0);
     90:   }
     91:   VecTaggerDestroy(&(*adaptor)->refineTag);
     92:   VecTaggerDestroy(&(*adaptor)->coarsenTag);
     93:   PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);
     94:   PetscHeaderDestroy(adaptor);
     95:   return(0);
     96: }

     98: /*@
     99:   DMAdaptorSetFromOptions - Sets a DMAdaptor object from options

    101:   Collective on DMAdaptor

    103:   Input Parameters:
    104: . adaptor - The DMAdaptor object

    106:   Options Database Keys:
    107: + -adaptor_monitor <bool>        : Monitor the adaptation process
    108: . -adaptor_sequence_num <num>    : Number of adaptations to generate an optimal grid
    109: . -adaptor_target_num <num>      : Set the target number of vertices N_adapt, -1 for automatic determination
    110: . -adaptor_refinement_factor <r> : Set r such that N_adapt = r^dim N_orig
    111: . -adaptor_metric_h_min <min>    : Set the minimum eigenvalue of Hessian (sqr max edge length)
    112: - -adaptor_metric_h_max <max>    : Set the maximum eigenvalue of Hessian (sqr min edge length)

    114:   Level: beginner

    116: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
    117: @*/
    118: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
    119: {

    123:   PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");
    124:   PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
    125:   PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
    126:   PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
    127:   PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
    128:   PetscOptionsReal("-adaptor_metric_h_min", "Set the minimum eigenvalue of Hessian (sqr max edge length)", "DMAdaptor", adaptor->h_min, &adaptor->h_min, NULL);
    129:   PetscOptionsReal("-adaptor_metric_h_max", "Set the maximum eigenvalue of Hessian (sqr min edge length)", "DMAdaptor", adaptor->h_max, &adaptor->h_max, NULL);
    130:   PetscOptionsEnd();
    131:   VecTaggerSetFromOptions(adaptor->refineTag);
    132:   VecTaggerSetFromOptions(adaptor->coarsenTag);
    133:   return(0);
    134: }

    136: /*@
    137:   DMAdaptorView - Views a DMAdaptor object

    139:   Collective on DMAdaptor

    141:   Input Parameters:
    142: + adaptor     - The DMAdaptor object
    143: - viewer - The PetscViewer object

    145:   Level: beginner

    147: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
    148: @*/
    149: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
    150: {

    154:   PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
    155:   PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
    156:   PetscViewerASCIIPrintf(viewer, "  sequence length: %D\n", adaptor->numSeq);
    157:   VecTaggerView(adaptor->refineTag,  viewer);
    158:   VecTaggerView(adaptor->coarsenTag, viewer);
    159:   return(0);
    160: }

    162: /*@
    163:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

    165:   Not collective

    167:   Input Parameter:
    168: . adaptor   - The DMAdaptor object

    170:   Output Parameter:
    171: . snes - The solver

    173:   Level: intermediate

    175: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
    176: @*/
    177: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
    178: {
    182:   *snes = adaptor->snes;
    183:   return(0);
    184: }

    186: /*@
    187:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

    189:   Not collective

    191:   Input Parameters:
    192: + adaptor   - The DMAdaptor object
    193: - snes - The solver

    195:   Level: intermediate

    197:   Note: The solver MUST have an attached DM/DS, so that we know the exact solution

    199: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
    200: @*/
    201: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
    202: {

    208:   adaptor->snes = snes;
    209:   SNESGetDM(adaptor->snes, &adaptor->idm);
    210:   return(0);
    211: }

    213: /*@
    214:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations

    216:   Not collective

    218:   Input Parameter:
    219: . adaptor - The DMAdaptor object

    221:   Output Parameter:
    222: . num - The number of adaptations

    224:   Level: intermediate

    226: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
    227: @*/
    228: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
    229: {
    233:   *num = adaptor->numSeq;
    234:   return(0);
    235: }

    237: /*@
    238:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

    240:   Not collective

    242:   Input Parameters:
    243: + adaptor - The DMAdaptor object
    244: - num - The number of adaptations

    246:   Level: intermediate

    248: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
    249: @*/
    250: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
    251: {
    254:   adaptor->numSeq = num;
    255:   return(0);
    256: }

    258: /*@
    259:   DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity

    261:   Collective on DMAdaptor

    263:   Input Parameters:
    264: . adaptor - The DMAdaptor object

    266:   Level: beginner

    268: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
    269: @*/
    270: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
    271: {
    272:   PetscDS        prob;
    273:   PetscInt       Nf, f;

    277:   DMGetDS(adaptor->idm, &prob);
    278:   VecTaggerSetUp(adaptor->refineTag);
    279:   VecTaggerSetUp(adaptor->coarsenTag);
    280:   PetscDSGetNumFields(prob, &Nf);
    281:   PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
    282:   for (f = 0; f < Nf; ++f) {
    283:     PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
    284:     /* TODO Have a flag that forces projection rather than using the exact solution */
    285:     if (adaptor->exactSol[0]) {DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);}
    286:   }
    287:   return(0);
    288: }

    290: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
    291: {
    293:   *tfunc = adaptor->ops->transfersolution;
    294:   return(0);
    295: }

    297: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
    298: {
    300:   adaptor->ops->transfersolution = tfunc;
    301:   return(0);
    302: }

    304: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
    305: {
    306:   DM             plex;
    307:   PetscDS        prob;
    308:   PetscObject    obj;
    309:   PetscClassId   id;
    310:   PetscBool      isForest;

    314:   DMConvert(adaptor->idm, DMPLEX, &plex);
    315:   DMGetDS(adaptor->idm, &prob);
    316:   PetscDSGetDiscretization(prob, 0, &obj);
    317:   PetscObjectGetClassId(obj, &id);
    318:   DMIsForest(adaptor->idm, &isForest);
    319:   if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
    320:     if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
    321: #if defined(PETSC_HAVE_PRAGMATIC)
    322:     else          {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
    323: #else
    324:     else          {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
    325: #endif
    326:   }
    327:   if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
    328:   else                       {adaptor->femType = PETSC_TRUE;}
    329:   if (adaptor->femType) {
    330:     /* Compute local solution bc */
    331:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
    332:   } else {
    333:     PetscFV      fvm = (PetscFV) obj;
    334:     PetscLimiter noneLimiter;
    335:     Vec          grad;

    337:     PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
    338:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
    339:     /* Use no limiting when reconstructing gradients for adaptivity */
    340:     PetscFVGetLimiter(fvm, &adaptor->limiter);
    341:     PetscObjectReference((PetscObject) adaptor->limiter);
    342:     PetscLimiterCreate(PetscObjectComm((PetscObject) fvm), &noneLimiter);
    343:     PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
    344:     PetscFVSetLimiter(fvm, noneLimiter);
    345:     /* Get FVM data */
    346:     DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
    347:     VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
    348:     VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
    349:     /* Compute local solution bc */
    350:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
    351:     /* Compute gradients */
    352:     DMCreateGlobalVector(adaptor->gradDM, &grad);
    353:     DMPlexReconstructGradientsFVM(plex, locX, grad);
    354:     DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
    355:     DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
    356:     DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
    357:     VecDestroy(&grad);
    358:     VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
    359:   }
    360:   DMDestroy(&plex);
    361:   return(0);
    362: }

    364: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
    365: {
    366:   PetscReal      time = 0.0;
    367:   Mat            interp;
    368:   void          *ctx;

    372:   DMGetApplicationContext(dm, &ctx);
    373:   if (adaptor->ops->transfersolution) {
    374:     (*adaptor->ops->transfersolution)(adaptor, dm, x, adm, ax, ctx);
    375:   } else {
    376:     switch (adaptor->adaptCriterion) {
    377:     case DM_ADAPTATION_LABEL:
    378:       DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
    379:       break;
    380:     case DM_ADAPTATION_REFINE:
    381:     case DM_ADAPTATION_METRIC:
    382:       DMCreateInterpolation(dm, adm, &interp, NULL);
    383:       MatInterpolate(interp, x, ax);
    384:       DMInterpolate(dm, interp, adm);
    385:       MatDestroy(&interp);
    386:       break;
    387:     default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %D", adaptor->adaptCriterion);
    388:     }
    389:   }
    390:   return(0);
    391: }

    393: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
    394: {
    395:   PetscDS        prob;
    396:   PetscObject    obj;
    397:   PetscClassId   id;

    401:   DMGetDS(adaptor->idm, &prob);
    402:   PetscDSGetDiscretization(prob, 0, &obj);
    403:   PetscObjectGetClassId(obj, &id);
    404:   if (id == PETSCFV_CLASSID) {
    405:     PetscFV fvm = (PetscFV) obj;

    407:     PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
    408:     /* Restore original limiter */
    409:     PetscFVSetLimiter(fvm, adaptor->limiter);

    411:     VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
    412:     VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
    413:     DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
    414:   }
    415:   return(0);
    416: }

    418: static PetscErrorCode DMAdaptorModifyHessian_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscScalar Hp[])
    419: {
    420:   PetscScalar   *Hpos;
    421:   PetscReal     *eigs;
    422:   PetscInt       i, j, k;

    426:   PetscMalloc2(dim*dim, &Hpos, dim, &eigs);
    427: #if 0
    428:   PetscPrintf(PETSC_COMM_SELF, "H = [");
    429:   for (i = 0; i < dim; ++i) {
    430:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "     ");}
    431:     for (j = 0; j < dim; ++j) {
    432:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    433:       PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
    434:     }
    435:     PetscPrintf(PETSC_COMM_SELF, "\n");
    436:   }
    437:   PetscPrintf(PETSC_COMM_SELF, "]\n");
    438: #endif
    439:   /* Symmetrize */
    440:   for (i = 0; i < dim; ++i) {
    441:     Hpos[i*dim+i] = Hp[i*dim+i];
    442:     for (j = i+1; j < dim; ++j) {
    443:       Hpos[i*dim+j] = 0.5*(Hp[i*dim+j] + Hp[j*dim+i]);
    444:       Hpos[j*dim+i] = Hpos[i*dim+j];
    445:     }
    446:   }
    447: #if 0
    448:   PetscPrintf(PETSC_COMM_SELF, "Hs = [");
    449:   for (i = 0; i < dim; ++i) {
    450:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "      ");}
    451:     for (j = 0; j < dim; ++j) {
    452:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    453:       PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
    454:     }
    455:     PetscPrintf(PETSC_COMM_SELF, "\n");
    456:   }
    457:   PetscPrintf(PETSC_COMM_SELF, "]\n");
    458: #endif
    459:   /* Compute eigendecomposition */
    460:   {
    461:     PetscScalar  *work;
    462:     PetscBLASInt lwork;

    464:     lwork = 5*dim;
    465:     PetscMalloc1(5*dim, &work);
    466:     {
    467:       PetscBLASInt lierr;
    468:       PetscBLASInt nb;

    470:       PetscBLASIntCast(dim, &nb);
    471:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    472: #if defined(PETSC_USE_COMPLEX)
    473:       {
    474:         PetscReal *rwork;
    475:         PetscMalloc1(3*dim, &rwork);
    476:         PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,rwork,&lierr));
    477:         PetscFree(rwork);
    478:       }
    479: #else
    480:       PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,&lierr));
    481: #endif
    482:       if (lierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
    483:       PetscFPTrapPop();
    484:     }
    485:     PetscFree(work);
    486:   }
    487: #if 0
    488:   PetscPrintf(PETSC_COMM_SELF, "L = [");
    489:   for (i = 0; i < dim; ++i) {
    490:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    491:     PetscPrintf(PETSC_COMM_SELF, "%g", eigs[i]);
    492:   }
    493:   PetscPrintf(PETSC_COMM_SELF, "]\n");
    494:   PetscPrintf(PETSC_COMM_SELF, "Q = [");
    495:   for (i = 0; i < dim; ++i) {
    496:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "     ");}
    497:     for (j = 0; j < dim; ++j) {
    498:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    499:       PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
    500:     }
    501:     PetscPrintf(PETSC_COMM_SELF, "\n");
    502:   }
    503:   PetscPrintf(PETSC_COMM_SELF, "]\n");
    504: #endif
    505:   /* Reflect to positive orthant, enforce maximum and minimum size, \lambda \propto 1/h^2
    506:        TODO get domain bounding box */
    507:   for (i = 0; i < dim; ++i) eigs[i] = PetscMin(h_max, PetscMax(h_min, PetscAbsReal(eigs[i])));
    508:   /* Reconstruct Hessian */
    509:   for (i = 0; i < dim; ++i) {
    510:     for (j = 0; j < dim; ++j) {
    511:       Hp[i*dim+j] = 0.0;
    512:       for (k = 0; k < dim; ++k) {
    513:         Hp[i*dim+j] += Hpos[k*dim+i] * eigs[k] * Hpos[k*dim+j];
    514:       }
    515:     }
    516:   }
    517:   PetscFree2(Hpos, eigs);
    518: #if 0
    519:   PetscPrintf(PETSC_COMM_SELF, "H+ = [");
    520:   for (i = 0; i < dim; ++i) {
    521:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "      ");}
    522:     for (j = 0; j < dim; ++j) {
    523:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
    524:       PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
    525:     }
    526:     PetscPrintf(PETSC_COMM_SELF, "\n");
    527:   }
    528:   PetscPrintf(PETSC_COMM_SELF, "]\n");
    529: #endif
    530:   return(0);
    531: }

    533: static void detHFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    534:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    535:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    536:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    537: {
    538:   const PetscInt p = 1;
    539:   PetscReal      detH = 0.0;

    541:   if      (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, u);
    542:   else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, u);
    543:   f0[0] = PetscPowReal(detH, p/(2.*p + dim));
    544: }

    546: /*
    547:   DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator

    549:   Input Parameters:
    550: + adaptor  - The DMAdaptor object
    551: . dim      - The topological dimension
    552: . cell     - The cell
    553: . field    - The field integrated over the cell
    554: . gradient - The gradient integrated over the cell
    555: . cg       - A PetscFVCellGeom struct
    556: - ctx      - A user context

    558:   Output Parameter:
    559: . errInd   - The error indicator

    561: .seealso: DMAdaptorComputeErrorIndicator()
    562: */
    563: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
    564: {
    565:   PetscReal err = 0.;
    566:   PetscInt  c, d;

    569:   for (c = 0; c < Nc; c++) {
    570:     for (d = 0; d < dim; ++d) {
    571:       err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
    572:     }
    573:   }
    574:   *errInd = cg->volume * err;
    575:   return(0);
    576: }

    578: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
    579: {
    580:   PetscDS         prob;
    581:   PetscObject     obj;
    582:   PetscClassId    id;
    583:   void           *ctx;
    584:   PetscQuadrature quad;
    585:   PetscInt        dim, d, cdim, Nc;
    586:   PetscErrorCode  ierr;

    589:   *errInd = 0.;
    590:   DMGetDimension(plex, &dim);
    591:   DMGetCoordinateDim(plex, &cdim);
    592:   DMGetApplicationContext(plex, &ctx);
    593:   DMGetDS(plex, &prob);
    594:   PetscDSGetDiscretization(prob, 0, &obj);
    595:   PetscObjectGetClassId(obj, &id);
    596:   if (id == PETSCFV_CLASSID) {
    597:     const PetscScalar *pointSols;
    598:     const PetscScalar *pointSol;
    599:     const PetscScalar *pointGrad;
    600:     PetscFVCellGeom   *cg;

    602:     PetscFVGetNumComponents((PetscFV) obj, &Nc);
    603:     VecGetArrayRead(locX, &pointSols);
    604:     DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
    605:     DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
    606:     DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
    607:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
    608:     VecRestoreArrayRead(locX, &pointSols);
    609:   } else {
    610:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
    611:     PetscFVCellGeom  cg;
    612:     PetscFEGeom      fegeom;
    613:     const PetscReal *quadWeights;
    614:     PetscReal       *coords;
    615:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

    617:     fegeom.dim      = dim;
    618:     fegeom.dimEmbed = cdim;
    619:     PetscDSGetNumFields(prob, &Nf);
    620:     PetscFEGetQuadrature((PetscFE) obj, &quad);
    621:     DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
    622:     PetscFEGetDimension((PetscFE) obj, &Nb);
    623:     PetscFEGetNumComponents((PetscFE) obj, &Nc);
    624:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
    625:     PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);
    626:     PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
    627:     DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
    628:     DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
    629:     PetscArrayzero(gradient, cdim*Nc);
    630:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
    631:       PetscInt qc = 0, q;

    633:       PetscDSGetDiscretization(prob, f, &obj);
    634:       PetscArrayzero(interpolant,Nc);
    635:       PetscArrayzero(interpolantGrad, cdim*Nc);
    636:       for (q = 0; q < Nq; ++q) {
    637:         PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, 1, x, &fegeom, q, interpolant, interpolantGrad);
    638:         for (fc = 0; fc < Nc; ++fc) {
    639:           const PetscReal wt = quadWeights[q*qNc+qc+fc];

    641:           field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
    642:           for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
    643:         }
    644:       }
    645:       fieldOffset += Nb;
    646:       qc          += Nc;
    647:     }
    648:     PetscFree2(interpolant, interpolantGrad);
    649:     DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
    650:     for (fc = 0; fc < Nc; ++fc) {
    651:       field[fc] /= cg.volume;
    652:       for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
    653:     }
    654:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);
    655:     PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
    656:   }
    657:   return(0);
    658: }

    660: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
    661: {
    662:   PetscDS        prob;
    663:   void          *ctx;
    664:   MPI_Comm       comm;
    665:   PetscInt       numAdapt = adaptor->numSeq, adaptIter;
    666:   PetscInt       dim, coordDim, numFields, cStart, cEnd, c;

    670:   DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
    671:   VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
    672:   PetscObjectGetComm((PetscObject) adaptor, &comm);
    673:   DMGetDimension(adaptor->idm, &dim);
    674:   DMGetCoordinateDim(adaptor->idm, &coordDim);
    675:   DMGetApplicationContext(adaptor->idm, &ctx);
    676:   DMGetDS(adaptor->idm, &prob);
    677:   PetscDSGetNumFields(prob, &numFields);

    679:   /* Adapt until nothing changes */
    680:   /* Adapt for a specified number of iterates */
    681:   for (adaptIter = 0; adaptIter < numAdapt-1; ++adaptIter) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));}
    682:   for (adaptIter = 0; adaptIter < numAdapt;   ++adaptIter) {
    683:     PetscBool adapted = PETSC_FALSE;
    684:     DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
    685:     Vec       x       = adaptIter ? *ax  : inx, locX, ox;

    687:     DMGetLocalVector(dm, &locX);
    688:     DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
    689:     DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
    690:     DMAdaptorPreAdapt(adaptor, locX);
    691:     if (doSolve) {
    692:       SNES snes;

    694:       DMAdaptorGetSolver(adaptor, &snes);
    695:       SNESSolve(snes, NULL, adaptIter ? *ax : x);
    696:     }
    697:     /* DMAdaptorMonitor(adaptor);
    698:        Print iterate, memory used, DM, solution */
    699:     switch (adaptor->adaptCriterion) {
    700:     case DM_ADAPTATION_REFINE:
    701:       DMRefine(dm, comm, &odm);
    702:       if (!odm) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
    703:       adapted = PETSC_TRUE;
    704:       break;
    705:     case DM_ADAPTATION_LABEL:
    706:     {
    707:       /* Adapt DM
    708:            Create local solution
    709:            Reconstruct gradients (FVM) or solve adjoint equation (FEM)
    710:            Produce cellwise error indicator */
    711:       DM                 plex;
    712:       DMLabel            adaptLabel;
    713:       IS                 refineIS, coarsenIS;
    714:       Vec                errVec;
    715:       PetscScalar       *errArray;
    716:       const PetscScalar *pointSols;
    717:       PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
    718:       PetscInt           nRefine, nCoarsen;

    720:       DMConvert(dm, DMPLEX, &plex);
    721:       DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
    722:       DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);

    724:       VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
    725:       VecSetUp(errVec);
    726:       VecGetArray(errVec, &errArray);
    727:       VecGetArrayRead(locX, &pointSols);
    728:       for (c = cStart; c < cEnd; ++c) {
    729:         PetscReal errInd;

    731:         DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);
    732:         errArray[c-cStart] = errInd;
    733:         minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
    734:         minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
    735:       }
    736:       VecRestoreArrayRead(locX, &pointSols);
    737:       VecRestoreArray(errVec, &errArray);
    738:       PetscGlobalMinMaxReal(PetscObjectComm((PetscObject) adaptor), minMaxInd, minMaxIndGlobal);
    739:       PetscInfo2(adaptor, "DMAdaptor: error indicator range (%E, %E)\n", minMaxIndGlobal[0], minMaxIndGlobal[1]);
    740:       /*     Compute IS from VecTagger */
    741:       VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS);
    742:       VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS);
    743:       ISGetSize(refineIS, &nRefine);
    744:       ISGetSize(coarsenIS, &nCoarsen);
    745:       PetscInfo2(adaptor, "DMAdaptor: numRefine %D, numCoarsen %D\n", nRefine, nCoarsen);
    746:       if (nRefine)  {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE,  refineIS);}
    747:       if (nCoarsen) {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);}
    748:       ISDestroy(&coarsenIS);
    749:       ISDestroy(&refineIS);
    750:       VecDestroy(&errVec);
    751:       /*     Adapt DM from label */
    752:       if (nRefine || nCoarsen) {
    753:         DMAdaptLabel(dm, adaptLabel, &odm);
    754:         adapted = PETSC_TRUE;
    755:       }
    756:       DMLabelDestroy(&adaptLabel);
    757:       DMDestroy(&plex);
    758:     }
    759:     break;
    760:     case DM_ADAPTATION_METRIC:
    761:     {
    762:       DM           dmGrad,   dmHess,   dmMetric;
    763:       PetscDS      probGrad, probHess;
    764:       Vec          xGrad,    xHess,    metric;
    765:       PetscSection sec, msec;
    766:       PetscScalar *H, *M, integral;
    767:       PetscReal    N;
    768:       DMLabel      bdLabel;
    769:       PetscInt     Nd = coordDim*coordDim, f, vStart, vEnd, v;

    771:       /*     Compute vertexwise gradients from cellwise gradients */
    772:       DMClone(dm, &dmGrad);
    773:       DMClone(dm, &dmHess);
    774:       DMGetDS(dmGrad, &probGrad);
    775:       DMGetDS(dmHess, &probHess);
    776:       for (f = 0; f < numFields; ++f) {
    777:         PetscFE         fe, feGrad, feHess;
    778:         PetscDualSpace  Q;
    779:         DM              K;
    780:         PetscQuadrature q;
    781:         PetscInt        Nc, qorder;
    782:         const char     *prefix;

    784:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
    785:         PetscFEGetNumComponents(fe, &Nc);
    786:         PetscFEGetDualSpace(fe, &Q);
    787:         PetscDualSpaceGetDM(Q, &K);
    788:         DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
    789:         PetscFEGetQuadrature(fe, &q);
    790:         PetscQuadratureGetOrder(q, &qorder);
    791:         PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);
    792:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feGrad);
    793:         PetscDSSetDiscretization(probGrad, f, (PetscObject) feGrad);
    794:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feHess);
    795:         PetscDSSetDiscretization(probHess, f, (PetscObject) feHess);
    796:         PetscFEDestroy(&feGrad);
    797:         PetscFEDestroy(&feHess);
    798:       }
    799:       DMGetGlobalVector(dmGrad, &xGrad);
    800:       VecViewFromOptions(x, NULL, "-sol_adapt_loc_pre_view");
    801:       DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
    802:       VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
    803:       /*     Compute vertexwise Hessians from cellwise Hessians */
    804:       DMGetGlobalVector(dmHess, &xHess);
    805:       DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
    806:       VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
    807:       /*     Compute metric */
    808:       DMClone(dm, &dmMetric);
    809:       DMGetLocalSection(dm, &sec);
    810:       DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    811:       PetscSectionCreate(PetscObjectComm((PetscObject) dm), &msec);
    812:       PetscSectionSetNumFields(msec, 1);
    813:       PetscSectionSetFieldComponents(msec, 0, Nd);
    814:       PetscSectionSetChart(msec, vStart, vEnd);
    815:       for (v = vStart; v < vEnd; ++v) {
    816:         PetscSectionSetDof(msec, v, Nd);
    817:         PetscSectionSetFieldDof(msec, v, 0, Nd);
    818:       }
    819:       PetscSectionSetUp(msec);
    820:       DMSetLocalSection(dmMetric, msec);
    821:       PetscSectionDestroy(&msec);
    822:       DMGetLocalVector(dmMetric, &metric);
    823:       /*       N is the target size */
    824:       N    = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
    825:       if (adaptor->monitor) {PetscPrintf(PETSC_COMM_SELF, "N_orig: %D N_adapt: %g\n", vEnd - vStart, N);}
    826:       /*       |H| means take the absolute value of eigenvalues */
    827:       VecGetArray(xHess, &H);
    828:       VecGetArray(metric, &M);
    829:       for (v = vStart; v < vEnd; ++v) {
    830:         PetscScalar *Hp;

    832:         DMPlexPointLocalRef(dmHess, v, H, &Hp);
    833:         DMAdaptorModifyHessian_Private(coordDim, adaptor->h_min, adaptor->h_max, Hp);
    834:       }
    835:       /*       Pointwise on vertices M(x) = N^{2/d} (\int_\Omega det(|H|)^{p/(2p+d)})^{-2/d} det(|H|)^{-1/(2p+d)} |H| for L_p */
    836:       PetscDSSetObjective(probHess, 0, detHFunc);
    837:       DMPlexComputeIntegralFEM(dmHess, xHess, &integral, NULL);
    838:       for (v = vStart; v < vEnd; ++v) {
    839:         const PetscInt     p = 1;
    840:         const PetscScalar *Hp;
    841:         PetscScalar       *Mp;
    842:         PetscReal          detH, fact;
    843:         PetscInt           i;

    845:         DMPlexPointLocalRead(dmHess, v, H, (void *) &Hp);
    846:         DMPlexPointLocalRef(dmMetric, v, M, &Mp);
    847:         if      (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, Hp);
    848:         else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, Hp);
    849:         else SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "Dimension %d not supported", dim);
    850:         fact = PetscPowReal(N, 2.0/dim) * PetscPowReal(PetscRealPart(integral), -2.0/dim) * PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim));
    851: #if 0
    852:         PetscPrintf(PETSC_COMM_SELF, "fact: %g integral: %g |detH|: %g termA: %g termB: %g\n", fact, integral, PetscAbsReal(detH), PetscPowReal(integral, -2.0/dim), PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim)));
    853:         DMPrintCellMatrix(v, "H", coordDim, coordDim, Hp);
    854: #endif
    855:         for (i = 0; i < Nd; ++i) {
    856:           Mp[i] = fact * Hp[i];
    857:         }
    858:       }
    859:       VecRestoreArray(xHess, &H);
    860:       VecRestoreArray(metric, &M);
    861:       /*     Adapt DM from metric */
    862:       DMGetLabel(dm, "marker", &bdLabel);
    863:       DMAdaptMetric(dm, metric, bdLabel, &odm);
    864:       adapted = PETSC_TRUE;
    865:       /* Cleanup */
    866:       DMRestoreLocalVector(dmMetric, &metric);
    867:       DMDestroy(&dmMetric);
    868:       DMRestoreGlobalVector(dmHess, &xHess);
    869:       DMRestoreGlobalVector(dmGrad, &xGrad);
    870:       DMDestroy(&dmGrad);
    871:       DMDestroy(&dmHess);
    872:     }
    873:     break;
    874:     default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
    875:     }
    876:     DMAdaptorPostAdapt(adaptor);
    877:     DMRestoreLocalVector(dm, &locX);
    878:     /* If DM was adapted, replace objects and recreate solution */
    879:     if (adapted) {
    880:       const char *name;

    882:       PetscObjectGetName((PetscObject) dm, &name);
    883:       PetscObjectSetName((PetscObject) odm, name);
    884:       /* Reconfigure solver */
    885:       SNESReset(adaptor->snes);
    886:       SNESSetDM(adaptor->snes, odm);
    887:       DMAdaptorSetSolver(adaptor, adaptor->snes);
    888:       DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
    889:       SNESSetFromOptions(adaptor->snes);
    890:       /* Transfer system */
    891:       DMCopyDisc(adaptor->idm, odm);
    892:       /* Transfer solution */
    893:       DMCreateGlobalVector(odm, &ox);
    894:       PetscObjectGetName((PetscObject) x, &name);
    895:       PetscObjectSetName((PetscObject) ox, name);
    896:       DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
    897:       /* Cleanup adaptivity info */
    898:       if (adaptIter > 0) {PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));}
    899:       DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
    900:       DMDestroy(&dm);
    901:       VecDestroy(&x);
    902:       *adm = odm;
    903:       *ax  = ox;
    904:     } else {
    905:       *adm = dm;
    906:       *ax  = x;
    907:       adaptIter = numAdapt;
    908:     }
    909:     if (adaptIter < numAdapt-1) {
    910:       DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
    911:       VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
    912:     }
    913:   }
    914:   DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
    915:   VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
    916:   return(0);
    917: }

    919: /*@
    920:   DMAdaptorAdapt - Creates a new DM that is adapted to the problem

    922:   Not collective

    924:   Input Parameter:
    925: + adaptor  - The DMAdaptor object
    926: . x        - The global approximate solution
    927: - strategy - The adaptation strategy

    929:   Output Parameters:
    930: + adm - The adapted DM
    931: - ax  - The adapted solution

    933:   Options database keys:
    934: + -snes_adapt <strategy> : initial, sequential, multigrid
    935: . -adapt_gradient_view : View the Clement interpolant of the solution gradient
    936: . -adapt_hessian_view : View the Clement interpolant of the solution Hessian
    937: - -adapt_metric_view : View the metric tensor for adaptive mesh refinement

    939:   Note: The available adaptation strategies are:
    940: $ 1) Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
    941: $ 2) Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
    942: $ 3) Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid

    944:   Level: intermediate

    946: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
    947: @*/
    948: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
    949: {

    953:   switch (strategy)
    954:   {
    955:   case DM_ADAPTATION_INITIAL:
    956:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
    957:     break;
    958:   case DM_ADAPTATION_SEQUENTIAL:
    959:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
    960:     break;
    961:   default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
    962:   }
    963:   return(0);
    964: }