Actual source code: tagger.c

      1: #include <petsc/private/vecimpl.h>

      3: /*@
      4:   VecTaggerCreate - create a `VecTagger` context.

      6:   Collective

      8:   Input Parameter:
      9: . comm - communicator on which the `VecTagger` will operate

     11:   Output Parameter:
     12: . tagger - new Vec tagger context

     14:   Level: advanced

     16:   Notes:
     17:   This object is used to control the tagging/selection of index sets based on the values in a
     18:   vector. This is used, for example, in adaptive simulations when aspects are selected for
     19:   refinement or coarsening. The primary intent is that the selected index sets are based purely
     20:   on the values in the vector, though implementations that do not follow this intent are
     21:   possible.

     23:   Once a `VecTagger` is created (`VecTaggerCreate()`), optionally modified by options
     24:   (`VecTaggerSetFromOptions()`), and set up (`VecTaggerSetUp()`), it is applied to vectors with
     25:   `VecTaggerComputeIS()` to compute the selected index sets.

     27:   Provided implementations support tagging based on a box/interval of values
     28:   (`VECTAGGERABSOLUTE`), based on a box of values of relative to the range of values present in
     29:   the vector (`VECTAGGERRELATIVE`), based on where values fall in the cumulative distribution
     30:   of values in the vector (`VECTAGGERCDF`), and based on unions (`VECTAGGEROR`) or
     31:   intersections (`VECTAGGERAND`) of other criteria.

     33: .seealso: `VecTagger`, `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
     34: @*/
     35: PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
     36: {
     37:   VecTagger b;

     39:   PetscFunctionBegin;
     40:   PetscAssertPointer(tagger, 2);
     41:   PetscCall(VecTaggerInitializePackage());

     43:   PetscCall(PetscHeaderCreate(b, VEC_TAGGER_CLASSID, "VecTagger", "Vec Tagger", "Vec", comm, VecTaggerDestroy, VecTaggerView));
     44:   b->blocksize   = 1;
     45:   b->invert      = PETSC_FALSE;
     46:   b->setupcalled = PETSC_FALSE;
     47:   *tagger        = b;
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: /*@
     52:   VecTaggerSetType - set the Vec tagger implementation

     54:   Collective

     56:   Input Parameters:
     57: + tagger - the `VecTagger` context
     58: - type   - a known method

     60:   Options Database Key:
     61: . -vec_tagger_type type - Sets the method; see `VecTaggerType`

     63:   Level: advanced

     65:   Notes:
     66:   See "include/petscvec.h" for available methods (for instance)
     67: +    `VECTAGGERABSOLUTE` - tag based on a box of values
     68: .    `VECTAGGERRELATIVE` - tag based on a box relative to the range of values present in the vector
     69: .    `VECTAGGERCDF`      - tag based on a box in the cumulative distribution of values present in the vector
     70: .    `VECTAGGEROR`       - tag based on the union of a set of `VecTagger` contexts
     71: -    `VECTAGGERAND`      - tag based on the intersection of a set of other `VecTagger` contexts

     73: .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
     74: @*/
     75: PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
     76: {
     77:   PetscBool match;
     78:   PetscErrorCode (*r)(VecTagger);

     80:   PetscFunctionBegin;
     82:   PetscAssertPointer(type, 2);

     84:   PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
     85:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

     87:   PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
     88:   PetscCheck(r, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested VecTagger type %s", type);
     89:   /* Destroy the previous private VecTagger context */
     90:   PetscTryTypeMethod(tagger, destroy);
     91:   PetscCall(PetscMemzero(tagger->ops, sizeof(*tagger->ops)));
     92:   PetscCall(PetscObjectChangeTypeName((PetscObject)tagger, type));
     93:   tagger->ops->create = r;
     94:   PetscCall((*r)(tagger));
     95:   PetscFunctionReturn(PETSC_SUCCESS);
     96: }

     98: /*@
     99:   VecTaggerGetType - Gets the `VecTaggerType` name (as a string) from the `VecTagger`.

    101:   Not Collective

    103:   Input Parameter:
    104: . tagger - The `VecTagger` context

    106:   Output Parameter:
    107: . type - The `VecTagger` type name

    109:   Level: advanced

    111: .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
    112: @*/
    113: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
    114: {
    115:   PetscFunctionBegin;
    117:   PetscAssertPointer(type, 2);
    118:   PetscCall(VecTaggerRegisterAll());
    119:   *type = ((PetscObject)tagger)->type_name;
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /*@
    124:   VecTaggerDestroy - destroy a `VecTagger` context

    126:   Collective

    128:   Input Parameter:
    129: . tagger - address of tagger

    131:   Level: advanced

    133: .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
    134: @*/
    135: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
    136: {
    137:   PetscFunctionBegin;
    138:   if (!*tagger) PetscFunctionReturn(PETSC_SUCCESS);
    140:   if (--((PetscObject)*tagger)->refct > 0) {
    141:     *tagger = NULL;
    142:     PetscFunctionReturn(PETSC_SUCCESS);
    143:   }
    144:   PetscTryTypeMethod(*tagger, destroy);
    145:   PetscCall(PetscHeaderDestroy(tagger));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*@
    150:   VecTaggerSetUp - set up a `VecTagger` context

    152:   Collective

    154:   Input Parameter:
    155: . tagger - Vec tagger object

    157:   Level: advanced

    159: .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`
    160: @*/
    161: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
    162: {
    163:   PetscFunctionBegin;
    164:   if (tagger->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
    165:   if (!((PetscObject)tagger)->type_name) PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
    166:   PetscTryTypeMethod(tagger, setup);
    167:   tagger->setupcalled = PETSC_TRUE;
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: /*@C
    172:   VecTaggerSetFromOptions - set `VecTagger` options using the options database

    174:   Logically Collective

    176:   Input Parameter:
    177: . tagger - vec tagger

    179:   Options Database Keys:
    180: + -vec_tagger_type       - implementation type, see `VecTaggerSetType()`
    181: . -vec_tagger_block_size - set the block size, see `VecTaggerSetBlockSize()`
    182: - -vec_tagger_invert     - invert the index set returned by `VecTaggerComputeIS()`

    184:   Level: advanced

    186: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
    187: @*/
    188: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
    189: {
    190:   VecTaggerType deft;
    191:   char          type[256];
    192:   PetscBool     flg;

    194:   PetscFunctionBegin;
    196:   PetscObjectOptionsBegin((PetscObject)tagger);
    197:   deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
    198:   PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
    199:   PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
    200:   PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
    201:   PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
    202:   PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
    203:   PetscOptionsEnd();
    204:   PetscFunctionReturn(PETSC_SUCCESS);
    205: }

    207: /*@
    208:   VecTaggerSetBlockSize - set the block size of the set of indices returned by `VecTaggerComputeIS()`.

    210:   Logically Collective

    212:   Input Parameters:
    213: + tagger    - vec tagger
    214: - blocksize - block size of the criteria used to tagger vectors

    216:   Level: advanced

    218:   Notes:
    219:   Values greater than one are useful when there are multiple criteria for determining which
    220:   indices to include in the set. For example, consider adaptive mesh refinement in a
    221:   multiphysics problem, with metrics of solution quality for multiple fields measure on each
    222:   cell. The size of the vector will be `[numCells` * numFields]`; the `VecTagger` block size
    223:   should be `numFields`; `VecTaggerComputeIS()` will return indices in the range `[0,
    224:   numCells)`, i.e., one index is given for each block of values.

    226:   Note that the block size of the vector does not have to match this block size.

    228: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
    229: @*/
    230: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
    231: {
    232:   PetscFunctionBegin;
    235:   tagger->blocksize = blocksize;
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /*@
    240:   VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.

    242:   Logically Collective

    244:   Input Parameter:
    245: . tagger - vec tagger

    247:   Output Parameter:
    248: . blocksize - block size of the vectors the tagger operates on

    250:   Level: advanced

    252: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
    253: @*/
    254: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
    255: {
    256:   PetscFunctionBegin;
    258:   PetscAssertPointer(blocksize, 2);
    259:   *blocksize = tagger->blocksize;
    260:   PetscFunctionReturn(PETSC_SUCCESS);
    261: }

    263: /*@
    264:   VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
    265:   then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
    266:   intersection of their exteriors.

    268:   Logically Collective

    270:   Input Parameters:
    271: + tagger - vec tagger
    272: - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is

    274:   Level: advanced

    276: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
    277: @*/
    278: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
    279: {
    280:   PetscFunctionBegin;
    283:   tagger->invert = invert;
    284:   PetscFunctionReturn(PETSC_SUCCESS);
    285: }

    287: /*@
    288:   VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted

    290:   Logically Collective

    292:   Input Parameter:
    293: . tagger - vec tagger

    295:   Output Parameter:
    296: . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is

    298:   Level: advanced

    300: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
    301: @*/
    302: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
    303: {
    304:   PetscFunctionBegin;
    306:   PetscAssertPointer(invert, 2);
    307:   *invert = tagger->invert;
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: /*@
    312:   VecTaggerView - view a `VecTagger` context

    314:   Collective

    316:   Input Parameters:
    317: + tagger - vec tagger
    318: - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`

    320:   Level: advanced

    322: .seealso: `VecTaggerCreate()`, `VecTagger`
    323: @*/
    324: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
    325: {
    326:   PetscBool isascii;

    328:   PetscFunctionBegin;
    330:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
    332:   PetscCheckSameComm(tagger, 1, viewer, 2);
    333:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    334:   if (isascii) {
    335:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
    336:     PetscCall(PetscViewerASCIIPushTab(viewer));
    337:     PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
    338:     PetscTryTypeMethod(tagger, view, viewer);
    339:     if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
    340:     PetscCall(PetscViewerASCIIPopTab(viewer));
    341:   }
    342:   PetscFunctionReturn(PETSC_SUCCESS);
    343: }

    345: /*@C
    346:   VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
    347:   in listed `PETSC_FALSE`

    349:   Collective

    351:   Input Parameters:
    352: + tagger - the `VecTagger` context
    353: - vec    - the vec to tag

    355:   Output Parameters:
    356: + numBoxes - the number of boxes in the tag definition
    357: . boxes    - a newly allocated list of boxes.  This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
    358: - listed   - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed

    360:   Level: advanced

    362:   Note:
    363:   A value is tagged if it is in any of the boxes, unless the tagger has been inverted (see `VecTaggerSetInvert()`/`VecTaggerGetInvert()`), in which case a value is tagged if it is in none of the boxes.

    365: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
    366: @*/
    367: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox *boxes[], PetscBool *listed)
    368: {
    369:   PetscInt vls, tbs;

    371:   PetscFunctionBegin;
    374:   PetscAssertPointer(numBoxes, 3);
    375:   PetscAssertPointer(boxes, 4);
    376:   PetscCall(VecGetLocalSize(vec, &vls));
    377:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
    378:   PetscCheck(vls % tbs == 0, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_INCOMP, "vec local size %" PetscInt_FMT " is not a multiple of tagger block size %" PetscInt_FMT, vls, tbs);
    379:   if (tagger->ops->computeboxes) {
    380:     *listed = PETSC_TRUE;
    381:     PetscUseTypeMethod(tagger, computeboxes, vec, numBoxes, boxes, listed);
    382:   } else *listed = PETSC_FALSE;
    383:   PetscFunctionReturn(PETSC_SUCCESS);
    384: }

    386: /*@C
    387:   VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values

    389:   Collective

    391:   Input Parameters:
    392: + tagger - the `VecTagger` context
    393: - vec    - the vec to tag

    395:   Output Parameters:
    396: + is     - a list of the indices tagged by the tagger, i.e., if the number of local indices will be n / bs, where n is `VecGetLocalSize()` and bs is `VecTaggerGetBlockSize()`.
    397: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed

    399:   Level: advanced

    401: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
    402: @*/
    403: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS is[], PetscBool *listed)
    404: {
    405:   PetscInt vls, tbs;

    407:   PetscFunctionBegin;
    410:   PetscAssertPointer(is, 3);
    411:   PetscCall(VecGetLocalSize(vec, &vls));
    412:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
    413:   PetscCheck(vls % tbs == 0, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_INCOMP, "vec local size %" PetscInt_FMT " is not a multiple of tagger block size %" PetscInt_FMT, vls, tbs);
    414:   if (tagger->ops->computeis) {
    415:     PetscUseTypeMethod(tagger, computeis, vec, is, listed);
    416:   } else if (listed) *listed = PETSC_FALSE;
    417:   PetscFunctionReturn(PETSC_SUCCESS);
    418: }

    420: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
    421: {
    422:   PetscInt           numBoxes;
    423:   VecTaggerBox      *boxes;
    424:   PetscInt           numTagged, offset;
    425:   PetscInt          *tagged;
    426:   PetscInt           bs, b, i, j, k, n;
    427:   PetscBool          invert;
    428:   const PetscScalar *vecArray;
    429:   PetscBool          boxlisted;

    431:   PetscFunctionBegin;
    432:   PetscCall(VecTaggerGetBlockSize(tagger, &bs));
    433:   PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
    434:   if (!boxlisted) {
    435:     if (listed) *listed = PETSC_FALSE;
    436:     PetscFunctionReturn(PETSC_SUCCESS);
    437:   }
    438:   PetscCall(VecGetArrayRead(vec, &vecArray));
    439:   PetscCall(VecGetLocalSize(vec, &n));
    440:   invert    = tagger->invert;
    441:   numTagged = 0;
    442:   offset    = 0;
    443:   tagged    = NULL;
    444:   PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
    445:   n /= bs;
    446:   for (i = 0; i < 2; i++) {
    447:     if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
    448:     for (j = 0; j < n; j++) {
    449:       for (k = 0; k < numBoxes; k++) {
    450:         for (b = 0; b < bs; b++) {
    451:           PetscScalar  val = vecArray[j * bs + b];
    452:           PetscInt     l   = k * bs + b;
    453:           VecTaggerBox box;
    454:           PetscBool    in;

    456:           box = boxes[l];
    457: #if !defined(PETSC_USE_COMPLEX)
    458:           in = (PetscBool)((box.min <= val) && (val <= box.max));
    459: #else
    460:           in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
    461: #endif
    462:           if (!in) break;
    463:         }
    464:         if (b == bs) break;
    465:       }
    466:       if ((PetscBool)(k < numBoxes) ^ invert) {
    467:         if (!i) numTagged++;
    468:         else tagged[offset++] = j;
    469:       }
    470:     }
    471:   }
    472:   PetscCall(VecRestoreArrayRead(vec, &vecArray));
    473:   PetscCall(PetscFree(boxes));
    474:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
    475:   PetscCall(ISSort(*is));
    476:   if (listed) *listed = PETSC_TRUE;
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }