Actual source code: tagger.c
1: #include <petsc/private/vecimpl.h>
3: /*@C
4: VecTaggerCreate - create a `VecTagger` context. This object is used to control the tagging/selection of index sets
5: based on the values in a vector. This is used, for example, in adaptive simulations when aspects are selected for
6: refinement or coarsening. The primary intent is that the selected index sets are based purely on the values in the
7: vector, though implementations that do not follow this intent are possible.
9: Once a `VecTagger` is created (`VecTaggerCreate()`), optionally modified by options (`VecTaggerSetFromOptions()`), and
10: set up (`VecTaggerSetUp()`), it is applied to vectors with `VecTaggerComputeIS()` to compute the selected index sets.
12: In many cases, the selection criteria for an index is whether the corresponding value falls within a collection of
13: boxes: for this common case, `VecTaggerCreateBoxes()` can also be used to determine those boxes.
15: Provided implementations support tagging based on a box/interval of values (`VECTAGGERABSOLUTE`), based on a box of
16: values of relative to the range of values present in the vector (`VECTAGGERRELATIVE`), based on where values fall in
17: the cumulative distribution of values in the vector (`VECTAGGERCDF`), and based on unions (`VECTAGGEROR`) or
18: intersections (`VECTAGGERAND`) of other criteria.
20: Collective
22: Input Parameter:
23: . comm - communicator on which the `VecTagger` will operate
25: Output Parameter:
26: . tagger - new Vec tagger context
28: Level: advanced
30: .seealso: `VecTagger`, `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
31: @*/
32: PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
33: {
34: VecTagger b;
36: PetscFunctionBegin;
38: PetscCall(VecTaggerInitializePackage());
40: PetscCall(PetscHeaderCreate(b, VEC_TAGGER_CLASSID, "VecTagger", "Vec Tagger", "Vec", comm, VecTaggerDestroy, VecTaggerView));
42: b->blocksize = 1;
43: b->invert = PETSC_FALSE;
44: b->setupcalled = PETSC_FALSE;
46: *tagger = b;
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: /*@C
51: VecTaggerSetType - set the Vec tagger implementation
53: Collective
55: Input Parameters:
56: + tagger - the `VecTagger` context
57: - type - a known method
59: Options Database Key:
60: . -vec_tagger_type <type> - Sets the method; use -help for a list
61: of available methods (for instance, absolute, relative, cdf, or, and)
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;
84: PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
85: if (match) PetscFunctionReturn(PETSC_SUCCESS);
87: PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
88: PetscCheck(r, PETSC_COMM_SELF, 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: /*@C
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;
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()`, `VecTaggerSetUp()`
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()`
188: @*/
189: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
190: {
191: VecTaggerType deft;
192: char type[256];
193: PetscBool flg;
195: PetscFunctionBegin;
197: PetscObjectOptionsBegin((PetscObject)tagger);
198: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
199: PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
200: PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
201: PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
202: PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
203: PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
204: PetscOptionsEnd();
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@C
209: VecTaggerSetBlockSize - block size of the set of indices returned by `VecTaggerComputeIS()`. Values greater than one
210: are useful when there are multiple criteria for determining which indices to include in the set. For example,
211: consider adaptive mesh refinement in a multiphysics problem, with metrics of solution quality for multiple fields
212: measure on each cell. The size of the vector will be [numCells * numFields]; the VecTagger block size should be
213: numFields; VecTaggerComputeIS() will return indices in the range [0,numCells), i.e., one index is given for each
214: block of values.
216: Note that the block size of the vector does not have to match.
218: Note also that the index set created in `VecTaggerComputeIS()` has block size: it is an index set over the list of
219: items that the vector refers to, not to the vector itself.
221: Logically Collective
223: Input Parameters:
224: + tagger - vec tagger
225: - blocksize - block size of the criteria used to tagger vectors
227: Level: advanced
229: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
230: @*/
231: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
232: {
233: PetscFunctionBegin;
236: tagger->blocksize = blocksize;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@C
241: VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.
243: Logically Collective
245: Input Parameter:
246: . tagger - vec tagger
248: Output Parameter:
249: . blocksize - block size of the vectors the tagger operates on
251: Level: advanced
253: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
254: @*/
255: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
256: {
257: PetscFunctionBegin;
260: *blocksize = tagger->blocksize;
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@C
265: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
266: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
267: intersection of their exteriors.
269: Logically Collective
271: Input Parameters:
272: + tagger - vec tagger
273: - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
275: Level: advanced
277: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
278: @*/
279: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
280: {
281: PetscFunctionBegin;
284: tagger->invert = invert;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@C
289: VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted
291: Logically Collective
293: Input Parameter:
294: . tagger - vec tagger
296: Output Parameter:
297: . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
299: Level: advanced
301: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
302: @*/
303: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
304: {
305: PetscFunctionBegin;
308: *invert = tagger->invert;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: VecTaggerView - view a `VecTagger` context
315: Collective
317: Input Parameters:
318: + tagger - vec tagger
319: - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`
321: Level: advanced
323: .seealso: `VecTaggerCreate()`, `VecTagger`, `VecTaggerCreate()`
324: @*/
325: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
326: {
327: PetscBool iascii;
329: PetscFunctionBegin;
331: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
333: PetscCheckSameComm(tagger, 1, viewer, 2);
334: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
335: if (iascii) {
336: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
337: PetscCall(PetscViewerASCIIPushTab(viewer));
338: PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
339: PetscTryTypeMethod(tagger, view, viewer);
340: if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
341: PetscCall(PetscViewerASCIIPopTab(viewer));
342: }
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@C
347: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
348: in listed `PETSC_FALSE`
350: Collective on VecTagger
352: Input Parameters:
353: + tagger - the `VecTagger` context
354: - vec - the vec to tag
356: Output Parameters:
357: + numBoxes - the number of boxes in the tag definition
358: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
359: - listed - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed
361: Level: advanced
363: Note:
364: 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.
366: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
367: @*/
368: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox **boxes, PetscBool *listed)
369: {
370: PetscInt vls, tbs;
372: PetscFunctionBegin;
377: PetscCall(VecGetLocalSize(vec, &vls));
378: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
379: 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);
380: if (tagger->ops->computeboxes) {
381: *listed = PETSC_TRUE;
382: PetscCall((*tagger->ops->computeboxes)(tagger, vec, numBoxes, boxes, listed));
383: } else *listed = PETSC_FALSE;
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@C
388: VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values
390: Collective
392: Input Parameters:
393: + tagger - the `VecTagger` context
394: - vec - the vec to tag
396: Output Parameters:
397: + 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()`.
398: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed
400: Level: advanced
402: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
403: @*/
404: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
405: {
406: PetscInt vls, tbs;
408: PetscFunctionBegin;
412: PetscCall(VecGetLocalSize(vec, &vls));
413: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
414: 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);
415: if (tagger->ops->computeis) {
416: PetscCall((*tagger->ops->computeis)(tagger, vec, is, listed));
417: } else if (listed) *listed = PETSC_FALSE;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
422: {
423: PetscInt numBoxes;
424: VecTaggerBox *boxes;
425: PetscInt numTagged, offset;
426: PetscInt *tagged;
427: PetscInt bs, b, i, j, k, n;
428: PetscBool invert;
429: const PetscScalar *vecArray;
430: PetscBool boxlisted;
432: PetscFunctionBegin;
433: PetscCall(VecTaggerGetBlockSize(tagger, &bs));
434: PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
435: if (!boxlisted) {
436: if (listed) *listed = PETSC_FALSE;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
439: PetscCall(VecGetArrayRead(vec, &vecArray));
440: PetscCall(VecGetLocalSize(vec, &n));
441: invert = tagger->invert;
442: numTagged = 0;
443: offset = 0;
444: tagged = NULL;
445: PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
446: n /= bs;
447: for (i = 0; i < 2; i++) {
448: if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
449: for (j = 0; j < n; j++) {
450: for (k = 0; k < numBoxes; k++) {
451: for (b = 0; b < bs; b++) {
452: PetscScalar val = vecArray[j * bs + b];
453: PetscInt l = k * bs + b;
454: VecTaggerBox box;
455: PetscBool in;
457: box = boxes[l];
458: #if !defined(PETSC_USE_COMPLEX)
459: in = (PetscBool)((box.min <= val) && (val <= box.max));
460: #else
461: in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
462: #endif
463: if (!in) break;
464: }
465: if (b == bs) break;
466: }
467: if ((PetscBool)(k < numBoxes) ^ invert) {
468: if (!i) numTagged++;
469: else tagged[offset++] = j;
470: }
471: }
472: }
473: PetscCall(VecRestoreArrayRead(vec, &vecArray));
474: PetscCall(PetscFree(boxes));
475: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
476: PetscCall(ISSort(*is));
477: if (listed) *listed = PETSC_TRUE;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }