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; use -help for a list
     62:    of available methods (for instance, absolute, relative, cdf, or, and)

     64:   Level: advanced

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

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

     81:   PetscFunctionBegin;
     83:   PetscAssertPointer(type, 2);

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

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

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

    102:   Not Collective

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

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

    110:   Level: advanced

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

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

    127:   Collective

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

    132:   Level: advanced

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

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

    153:   Collective

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

    158:   Level: advanced

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

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

    175:   Logically Collective

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

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

    185:   Level: advanced

    187: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`

    189: @*/
    190: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
    191: {
    192:   VecTaggerType deft;
    193:   char          type[256];
    194:   PetscBool     flg;

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

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

    212:   Logically Collective

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

    218:   Level: advanced

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

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

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

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

    244:   Logically Collective

    246:   Input Parameter:
    247: . tagger - vec tagger

    249:   Output Parameter:
    250: . blocksize - block size of the vectors the tagger operates on

    252:   Level: advanced

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

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

    270:   Logically Collective

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

    276:   Level: advanced

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

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

    292:   Logically Collective

    294:   Input Parameter:
    295: . tagger - vec tagger

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

    300:   Level: advanced

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

    313: /*@
    314:   VecTaggerView - view a `VecTagger` context

    316:   Collective

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

    322:   Level: advanced

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

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

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

    351:   Collective

    353:   Input Parameters:
    354: + tagger - the `VecTagger` context
    355: - vec    - the vec to tag

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

    362:   Level: advanced

    364:   Note:
    365:   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.

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

    373:   PetscFunctionBegin;
    376:   PetscAssertPointer(numBoxes, 3);
    377:   PetscAssertPointer(boxes, 4);
    378:   PetscCall(VecGetLocalSize(vec, &vls));
    379:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
    380:   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);
    381:   if (tagger->ops->computeboxes) {
    382:     *listed = PETSC_TRUE;
    383:     PetscUseTypeMethod(tagger, computeboxes, vec, numBoxes, boxes, listed);
    384:   } else *listed = PETSC_FALSE;
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }

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

    391:   Collective

    393:   Input Parameters:
    394: + tagger - the `VecTagger` context
    395: - vec    - the vec to tag

    397:   Output Parameters:
    398: + 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()`.
    399: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed

    401:   Level: advanced

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

    409:   PetscFunctionBegin;
    412:   PetscAssertPointer(is, 3);
    413:   PetscCall(VecGetLocalSize(vec, &vls));
    414:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
    415:   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);
    416:   if (tagger->ops->computeis) {
    417:     PetscUseTypeMethod(tagger, computeis, vec, is, listed);
    418:   } else if (listed) *listed = PETSC_FALSE;
    419:   PetscFunctionReturn(PETSC_SUCCESS);
    420: }

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

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

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