1 #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
2
3 /*@
4 VecTaggerCreate - create a `VecTagger` context.
5
6 Collective
7
8 Input Parameter:
9 . comm - communicator on which the `VecTagger` will operate
10
11 Output Parameter:
12 . tagger - new Vec tagger context
13
14 Level: advanced
15
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.
22
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.
26
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.
32
33 .seealso: `VecTagger`, `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
34 @*/
VecTaggerCreate(MPI_Comm comm,VecTagger * tagger)35 PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
36 {
37 VecTagger b;
38
39 PetscFunctionBegin;
40 PetscAssertPointer(tagger, 2);
41 PetscCall(VecTaggerInitializePackage());
42
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 }
50
51 /*@
52 VecTaggerSetType - set the Vec tagger implementation
53
54 Collective
55
56 Input Parameters:
57 + tagger - the `VecTagger` context
58 - type - a known method
59
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)
63
64 Level: advanced
65
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
73
74 .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
75 @*/
VecTaggerSetType(VecTagger tagger,VecTaggerType type)76 PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
77 {
78 PetscBool match;
79 PetscErrorCode (*r)(VecTagger);
80
81 PetscFunctionBegin;
82 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
83 PetscAssertPointer(type, 2);
84
85 PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
86 if (match) PetscFunctionReturn(PETSC_SUCCESS);
87
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 }
98
99 /*@
100 VecTaggerGetType - Gets the `VecTaggerType` name (as a string) from the `VecTagger`.
101
102 Not Collective
103
104 Input Parameter:
105 . tagger - The `VecTagger` context
106
107 Output Parameter:
108 . type - The `VecTagger` type name
109
110 Level: advanced
111
112 .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
113 @*/
VecTaggerGetType(VecTagger tagger,VecTaggerType * type)114 PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
115 {
116 PetscFunctionBegin;
117 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
118 PetscAssertPointer(type, 2);
119 PetscCall(VecTaggerRegisterAll());
120 *type = ((PetscObject)tagger)->type_name;
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
124 /*@
125 VecTaggerDestroy - destroy a `VecTagger` context
126
127 Collective
128
129 Input Parameter:
130 . tagger - address of tagger
131
132 Level: advanced
133
134 .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
135 @*/
VecTaggerDestroy(VecTagger * tagger)136 PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
137 {
138 PetscFunctionBegin;
139 if (!*tagger) PetscFunctionReturn(PETSC_SUCCESS);
140 PetscValidHeaderSpecific(*tagger, VEC_TAGGER_CLASSID, 1);
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 }
149
150 /*@
151 VecTaggerSetUp - set up a `VecTagger` context
152
153 Collective
154
155 Input Parameter:
156 . tagger - Vec tagger object
157
158 Level: advanced
159
160 .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`
161 @*/
VecTaggerSetUp(VecTagger tagger)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 }
171
172 /*@C
173 VecTaggerSetFromOptions - set `VecTagger` options using the options database
174
175 Logically Collective
176
177 Input Parameter:
178 . tagger - vec tagger
179
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()`
184
185 Level: advanced
186
187 .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
188
189 @*/
VecTaggerSetFromOptions(VecTagger tagger)190 PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
191 {
192 VecTaggerType deft;
193 char type[256];
194 PetscBool flg;
195
196 PetscFunctionBegin;
197 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
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 }
208
209 /*@
210 VecTaggerSetBlockSize - set the block size of the set of indices returned by `VecTaggerComputeIS()`.
211
212 Logically Collective
213
214 Input Parameters:
215 + tagger - vec tagger
216 - blocksize - block size of the criteria used to tagger vectors
217
218 Level: advanced
219
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.
227
228 Note that the block size of the vector does not have to match this block size.
229
230 .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
231 @*/
VecTaggerSetBlockSize(VecTagger tagger,PetscInt blocksize)232 PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
233 {
234 PetscFunctionBegin;
235 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
236 PetscValidLogicalCollectiveInt(tagger, blocksize, 2);
237 tagger->blocksize = blocksize;
238 PetscFunctionReturn(PETSC_SUCCESS);
239 }
240
241 /*@
242 VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.
243
244 Logically Collective
245
246 Input Parameter:
247 . tagger - vec tagger
248
249 Output Parameter:
250 . blocksize - block size of the vectors the tagger operates on
251
252 Level: advanced
253
254 .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
255 @*/
VecTaggerGetBlockSize(VecTagger tagger,PetscInt * blocksize)256 PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
257 {
258 PetscFunctionBegin;
259 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
260 PetscAssertPointer(blocksize, 2);
261 *blocksize = tagger->blocksize;
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
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.
269
270 Logically Collective
271
272 Input Parameters:
273 + tagger - vec tagger
274 - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
275
276 Level: advanced
277
278 .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
279 @*/
VecTaggerSetInvert(VecTagger tagger,PetscBool invert)280 PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
281 {
282 PetscFunctionBegin;
283 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
284 PetscValidLogicalCollectiveBool(tagger, invert, 2);
285 tagger->invert = invert;
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
289 /*@
290 VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted
291
292 Logically Collective
293
294 Input Parameter:
295 . tagger - vec tagger
296
297 Output Parameter:
298 . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
299
300 Level: advanced
301
302 .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
303 @*/
VecTaggerGetInvert(VecTagger tagger,PetscBool * invert)304 PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
305 {
306 PetscFunctionBegin;
307 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
308 PetscAssertPointer(invert, 2);
309 *invert = tagger->invert;
310 PetscFunctionReturn(PETSC_SUCCESS);
311 }
312
313 /*@
314 VecTaggerView - view a `VecTagger` context
315
316 Collective
317
318 Input Parameters:
319 + tagger - vec tagger
320 - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`
321
322 Level: advanced
323
324 .seealso: `VecTaggerCreate()`, `VecTagger`
325 @*/
VecTaggerView(VecTagger tagger,PetscViewer viewer)326 PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
327 {
328 PetscBool isascii;
329
330 PetscFunctionBegin;
331 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
332 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
333 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
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 }
346
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`
350
351 Collective
352
353 Input Parameters:
354 + tagger - the `VecTagger` context
355 - vec - the vec to tag
356
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
361
362 Level: advanced
363
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.
366
367 .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
368 @*/
VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt * numBoxes,VecTaggerBox * boxes[],PetscBool * listed)369 PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox *boxes[], PetscBool *listed)
370 {
371 PetscInt vls, tbs;
372
373 PetscFunctionBegin;
374 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
375 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
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 }
387
388 /*@C
389 VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values
390
391 Collective
392
393 Input Parameters:
394 + tagger - the `VecTagger` context
395 - vec - the vec to tag
396
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
400
401 Level: advanced
402
403 .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
404 @*/
VecTaggerComputeIS(VecTagger tagger,Vec vec,IS is[],PetscBool * listed)405 PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS is[], PetscBool *listed)
406 {
407 PetscInt vls, tbs;
408
409 PetscFunctionBegin;
410 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID, 1);
411 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
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 }
421
VecTaggerComputeIS_FromBoxes(VecTagger tagger,Vec vec,IS * is,PetscBool * listed)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;
432
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;
457
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 }
481