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 iascii;

330:   PetscFunctionBegin;
332:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
334:   PetscCheckSameComm(tagger, 1, viewer, 2);
335:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
336:   if (iascii) {
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: }