1 #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
2
3 /*@C
4 VecTaggerCreate - create a Vec tagger 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.
8
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 comute the selected index sets.
11
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.
14
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.
19
20 Collective
21
22 Input Arguments:
23 . comm - communicator on which the vec tagger will operate
24
25 Output Arguments:
26 . tagger - new Vec tagger context
27
28 Level: advanced
29
30 .seealso: VecTaggerSetBlockSize(), VecTaggerSetFromOptions(), VecTaggerSetUp(), VecTaggerComputeIS(), VecTaggerComputeBoxes(), VecTaggerDestroy()
31 @*/
VecTaggerCreate(MPI_Comm comm,VecTagger * tagger)32 PetscErrorCode VecTaggerCreate(MPI_Comm comm,VecTagger *tagger)
33 {
34 PetscErrorCode ierr;
35 VecTagger b;
36
37 PetscFunctionBegin;
38 PetscValidPointer(tagger,2);
39 ierr = VecTaggerInitializePackage();CHKERRQ(ierr);
40
41 ierr = PetscHeaderCreate(b,VEC_TAGGER_CLASSID,"VecTagger","Vec Tagger","Vec",comm,VecTaggerDestroy,VecTaggerView);CHKERRQ(ierr);
42
43 b->blocksize = 1;
44 b->invert = PETSC_FALSE;
45 b->setupcalled = PETSC_FALSE;
46
47 *tagger = b;
48 PetscFunctionReturn(0);
49 }
50
51 /*@C
52 VecTaggerSetType - set the Vec tagger implementation
53
54 Collective on VecTagger
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 Notes:
65 See "include/petscvec.h" for available methods (for instance)
66 + VECTAGGERABSOLUTE - tag based on a box of values
67 . VECTAGGERRELATIVE - tag based on a box relative to the range of values present in the vector
68 . VECTAGGERCDF - tag based on a box in the cumulative distribution of values present in the vector
69 . VECTAGGEROR - tag based on the union of a set of VecTagger contexts
70 - VECTAGGERAND - tag based on the intersection of a set of other VecTagger contexts
71
72 Level: advanced
73
74 .seealso: VecTaggerType, VecTaggerCreate()
75 @*/
VecTaggerSetType(VecTagger tagger,VecTaggerType type)76 PetscErrorCode VecTaggerSetType(VecTagger tagger,VecTaggerType type)
77 {
78 PetscErrorCode ierr,(*r)(VecTagger);
79 PetscBool match;
80
81 PetscFunctionBegin;
82 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
83 PetscValidCharPointer(type,2);
84
85 ierr = PetscObjectTypeCompare((PetscObject)tagger,type,&match);CHKERRQ(ierr);
86 if (match) PetscFunctionReturn(0);
87
88 ierr = PetscFunctionListFind(VecTaggerList,type,&r);CHKERRQ(ierr);
89 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested VecTagger type %s",type);
90 /* Destroy the previous private VecTagger context */
91 if (tagger->ops->destroy) {
92 ierr = (*(tagger)->ops->destroy)(tagger);CHKERRQ(ierr);
93 }
94 ierr = PetscMemzero(tagger->ops,sizeof(*tagger->ops));CHKERRQ(ierr);
95 ierr = PetscObjectChangeTypeName((PetscObject)tagger,type);CHKERRQ(ierr);
96 tagger->ops->create = r;
97 ierr = (*r)(tagger);CHKERRQ(ierr);
98 PetscFunctionReturn(0);
99 }
100
101 /*@C
102 VecTaggerGetType - Gets the VecTagger type name (as a string) from the VecTagger.
103
104 Not Collective
105
106 Input Parameter:
107 . tagger - The Vec tagger context
108
109 Output Parameter:
110 . type - The VecTagger type name
111
112 Level: advanced
113
114 .seealso: VecTaggerSetType(), VecTaggerCreate()
115 @*/
VecTaggerGetType(VecTagger tagger,VecTaggerType * type)116 PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
117 {
118 PetscErrorCode ierr;
119
120 PetscFunctionBegin;
121 PetscValidHeaderSpecific(tagger, VEC_TAGGER_CLASSID,1);
122 PetscValidPointer(type,2);
123 ierr = VecTaggerRegisterAll();CHKERRQ(ierr);
124 *type = ((PetscObject)tagger)->type_name;
125 PetscFunctionReturn(0);
126 }
127
128 /*@
129 VecTaggerDestroy - destroy a VecTagger context
130
131 Collective
132
133 Input Arguments:
134 . tagger - address of tagger
135
136 Level: advanced
137
138 .seealso: VecTaggerCreate()
139 @*/
VecTaggerDestroy(VecTagger * tagger)140 PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
141 {
142 PetscErrorCode ierr;
143
144 PetscFunctionBegin;
145 if (!*tagger) PetscFunctionReturn(0);
146 PetscValidHeaderSpecific((*tagger),VEC_TAGGER_CLASSID,1);
147 if (--((PetscObject)(*tagger))->refct > 0) {*tagger = NULL; PetscFunctionReturn(0);}
148 if ((*tagger)->ops->destroy) {ierr = (*(*tagger)->ops->destroy)(*tagger);CHKERRQ(ierr);}
149 ierr = PetscHeaderDestroy(tagger);CHKERRQ(ierr);
150 PetscFunctionReturn(0);
151 }
152
153 /*@
154 VecTaggerSetUp - set up a VecTagger context
155
156 Collective
157
158 Input Arguments:
159 . tagger - Vec tagger object
160
161 Level: advanced
162
163 .seealso: VecTaggerSetFromOptions(), VecTaggerSetType()
164 @*/
VecTaggerSetUp(VecTagger tagger)165 PetscErrorCode VecTaggerSetUp(VecTagger tagger)
166 {
167 PetscErrorCode ierr;
168
169 PetscFunctionBegin;
170 if (tagger->setupcalled) PetscFunctionReturn(0);
171 if (!((PetscObject)tagger)->type_name) {ierr = VecTaggerSetType(tagger,VECTAGGERABSOLUTE);CHKERRQ(ierr);}
172 if (tagger->ops->setup) {ierr = (*tagger->ops->setup)(tagger);CHKERRQ(ierr);}
173 tagger->setupcalled = PETSC_TRUE;
174 PetscFunctionReturn(0);
175 }
176
177 /*@C
178 VecTaggerSetFromOptions - set VecTagger options using the options database
179
180 Logically Collective
181
182 Input Arguments:
183 . tagger - vec tagger
184
185 Options Database Keys:
186 + -vec_tagger_type - implementation type, see VecTaggerSetType()
187 . -vec_tagger_block_size - set the block size, see VecTaggerSetBlockSize()
188 - -vec_tagger_invert - invert the index set returned by VecTaggerComputeIS()
189
190 Level: advanced
191
192 @*/
VecTaggerSetFromOptions(VecTagger tagger)193 PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
194 {
195 VecTaggerType deft;
196 char type[256];
197 PetscErrorCode ierr;
198 PetscBool flg;
199
200 PetscFunctionBegin;
201 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
202 ierr = PetscObjectOptionsBegin((PetscObject)tagger);CHKERRQ(ierr);
203 deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
204 ierr = PetscOptionsFList("-vec_tagger_type","VecTagger implementation type","VecTaggerSetType",VecTaggerList,deft,type,256,&flg);CHKERRQ(ierr);
205 ierr = VecTaggerSetType(tagger,flg ? type : deft);CHKERRQ(ierr);
206 ierr = PetscOptionsInt("-vec_tagger_block_size","block size of the vectors the tagger operates on","VecTaggerSetBlockSize",tagger->blocksize,&tagger->blocksize,NULL);CHKERRQ(ierr);
207 ierr = PetscOptionsBool("-vec_tagger_invert","invert the set of indices returned by VecTaggerComputeIS()","VecTaggerSetInvert",tagger->invert,&tagger->invert,NULL);CHKERRQ(ierr);
208 if (tagger->ops->setfromoptions) {ierr = (*tagger->ops->setfromoptions)(PetscOptionsObject,tagger);CHKERRQ(ierr);}
209 ierr = PetscOptionsEnd();CHKERRQ(ierr);
210 PetscFunctionReturn(0);
211 }
212
213 /*@C
214 VecTaggerSetBlockSize - block size of the set of indices returned by VecTaggerComputeIS(). Values greater than one
215 are useful when there are multiple criteria for determining which indices to include in the set. For example,
216 consider adaptive mesh refinement in a multiphysics problem, with metrics of solution quality for multiple fields
217 measure on each cell. The size of the vector will be [numCells * numFields]; the VecTagger block size should be
218 numFields; VecTaggerComputeIS() will return indices in the range [0,numCells), i.e., one index is given for each
219 block of values.
220
221 Note that the block size of the vector does not have to match.
222
223 Note also that the index set created in VecTaggerComputeIS() has block size: it is an index set over the list of
224 items that the vector refers to, not to the vector itself.
225
226 Logically Collective
227
228 Input Arguments:
229 + tagger - vec tagger
230 - blocksize - block size of the criteria used to tagger vectors
231
232 Level: advanced
233
234 .seealso: VecTaggerComputeIS(), VecTaggerGetBlockSize(), VecSetBlockSize(), VecGetBlockSize()
235 @*/
VecTaggerSetBlockSize(VecTagger tagger,PetscInt blocksize)236 PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
237 {
238
239 PetscFunctionBegin;
240 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
241 PetscValidLogicalCollectiveInt(tagger,blocksize,2);
242 tagger->blocksize = blocksize;
243 PetscFunctionReturn(0);
244 }
245
246 /*@C
247 VecTaggerGetBlockSize - get the block size of the indices created by VecTaggerComputeIS().
248
249 Logically Collective
250
251 Input Arguments:
252 + tagger - vec tagger
253 - blocksize - block size of the vectors the tagger operates on
254
255 Level: advanced
256
257 .seealso: VecTaggerComputeIS(), VecTaggerSetBlockSize(),
258 @*/
VecTaggerGetBlockSize(VecTagger tagger,PetscInt * blocksize)259 PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
260 {
261
262 PetscFunctionBegin;
263 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
264 PetscValidPointer(blocksize,2);
265 *blocksize = tagger->blocksize;
266 PetscFunctionReturn(0);
267 }
268
269 /*@C
270 VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by VecTaggerComputeBoxes(),
271 then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
272 intersection of their exteriors.
273
274 Logically Collective
275
276 Input Arguments:
277 + tagger - vec tagger
278 - invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
279
280 Level: advanced
281
282 .seealso: VecTaggerComputeIS(), VecTaggerGetInvert()
283 @*/
VecTaggerSetInvert(VecTagger tagger,PetscBool invert)284 PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
285 {
286
287 PetscFunctionBegin;
288 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
289 PetscValidLogicalCollectiveBool(tagger,invert,2);
290 tagger->invert = invert;
291 PetscFunctionReturn(0);
292 }
293
294 /*@C
295 VecTaggerGetInvert - get whether the set of indices returned by VecTaggerComputeIS() are inverted
296
297 Logically Collective
298
299 Input Arguments:
300 + tagger - vec tagger
301 - invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
302
303 Level: advanced
304
305 .seealso: VecTaggerComputeIS(), VecTaggerSetInvert()
306 @*/
VecTaggerGetInvert(VecTagger tagger,PetscBool * invert)307 PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
308 {
309
310 PetscFunctionBegin;
311 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
312 PetscValidPointer(invert,2);
313 *invert = tagger->invert;
314 PetscFunctionReturn(0);
315 }
316
317 /*@C
318 VecTaggerView - view a VecTagger context
319
320 Collective
321
322 Input Arguments:
323 + tagger - vec tagger
324 - viewer - viewer to display tagger, for example PETSC_VIEWER_STDOUT_WORLD
325
326 Level: advanced
327
328 .seealso: VecTaggerCreate()
329 @*/
VecTaggerView(VecTagger tagger,PetscViewer viewer)330 PetscErrorCode VecTaggerView(VecTagger tagger,PetscViewer viewer)
331 {
332 PetscErrorCode ierr;
333 PetscBool iascii;
334
335 PetscFunctionBegin;
336 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
337 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger),&viewer);CHKERRQ(ierr);}
338 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
339 PetscCheckSameComm(tagger,1,viewer,2);
340 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
341 if (iascii) {
342 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)tagger,viewer);CHKERRQ(ierr);
343 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
344 ierr = PetscViewerASCIIPrintf(viewer,"Block size: %D\n",tagger->blocksize);CHKERRQ(ierr);
345 if (tagger->ops->view) {ierr = (*tagger->ops->view)(tagger,viewer);CHKERRQ(ierr);}
346 if (tagger->invert) {ierr = PetscViewerASCIIPrintf(viewer,"Inverting ISs.\n");CHKERRQ(ierr);}
347 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
348 }
349 PetscFunctionReturn(0);
350 }
351
352 /*@C
353 VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list
354
355 Collective on VecTagger
356
357 Input Aguments:
358 + tagger - the VecTagger context
359 - vec - the vec to tag
360
361 Output Arguments:
362 + numBoxes - the number of boxes in the tag definition
363 - boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * numBoxes) pairs that the user can free with PetscFree().
364
365 Notes:
366 . 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
368 Level: advanced
369
370 .seealso: VecTaggerComputeIS()
371 @*/
VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt * numBoxes,VecTaggerBox ** boxes)372 PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt *numBoxes,VecTaggerBox **boxes)
373 {
374 PetscInt vls, tbs;
375 PetscErrorCode ierr;
376
377 PetscFunctionBegin;
378 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
379 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
380 PetscValidIntPointer(numBoxes,3);
381 PetscValidPointer(boxes,4);
382 ierr = VecGetLocalSize(vec,&vls);CHKERRQ(ierr);
383 ierr = VecTaggerGetBlockSize(tagger,&tbs);CHKERRQ(ierr);
384 if (vls % tbs) SETERRQ2(PetscObjectComm((PetscObject)tagger),PETSC_ERR_ARG_INCOMP,"vec local size %D is not a multiple of tagger block size %D",vls,tbs);
385 if (tagger->ops->computeboxes) {ierr = (*tagger->ops->computeboxes) (tagger,vec,numBoxes,boxes);CHKERRQ(ierr);}
386 else {
387 const char *type;
388 ierr = PetscObjectGetType ((PetscObject)tagger,&type);CHKERRQ(ierr);
389 SETERRQ1(PetscObjectComm((PetscObject)tagger),PETSC_ERR_SUP,"VecTagger type %s does not compute value boxes",type);
390 }
391 PetscFunctionReturn(0);
392 }
393
394 /*@C
395 VecTaggerComputeIS - Use a VecTagger context to tag a set of indices based on a vector's values
396
397 Collective on VecTagger
398
399 Input Aguments:
400 + tagger - the VecTagger context
401 - vec - the vec to tag
402
403 Output Arguments:
404 . 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().
405
406 Level: advanced
407
408 .seealso: VecTaggerComputeBoxes()
409 @*/
VecTaggerComputeIS(VecTagger tagger,Vec vec,IS * is)410 PetscErrorCode VecTaggerComputeIS(VecTagger tagger,Vec vec,IS *is)
411 {
412 PetscInt vls, tbs;
413 PetscErrorCode ierr;
414
415 PetscFunctionBegin;
416 PetscValidHeaderSpecific(tagger,VEC_TAGGER_CLASSID,1);
417 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
418 PetscValidPointer(is,3);
419 ierr = VecGetLocalSize(vec,&vls);CHKERRQ(ierr);
420 ierr = VecTaggerGetBlockSize(tagger,&tbs);CHKERRQ(ierr);
421 if (vls % tbs) SETERRQ2(PetscObjectComm((PetscObject)tagger),PETSC_ERR_ARG_INCOMP,"vec local size %D is not a multiple of tagger block size %D",vls,tbs);
422 if (tagger->ops->computeis) {ierr = (*tagger->ops->computeis) (tagger,vec,is);CHKERRQ(ierr);}
423 else {
424 SETERRQ(PetscObjectComm((PetscObject)tagger),PETSC_ERR_SUP,"VecTagger type does not compute ISs");
425 }
426 PetscFunctionReturn(0);
427 }
428
VecTaggerComputeIS_FromBoxes(VecTagger tagger,Vec vec,IS * is)429 PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is)
430 { PetscInt numBoxes;
431 VecTaggerBox *boxes;
432 PetscInt numTagged, offset;
433 PetscInt *tagged;
434 PetscInt bs, b, i, j, k, n;
435 PetscBool invert;
436 const PetscScalar *vecArray;
437 PetscErrorCode ierr;
438
439 PetscFunctionBegin;
440 ierr = VecTaggerGetBlockSize(tagger,&bs);CHKERRQ(ierr);
441 ierr = VecTaggerComputeBoxes(tagger,vec,&numBoxes,&boxes);CHKERRQ(ierr);
442 ierr = VecGetArrayRead(vec, &vecArray);CHKERRQ(ierr);
443 ierr = VecGetLocalSize(vec, &n);CHKERRQ(ierr);
444 invert = tagger->invert;
445 numTagged = 0;
446 offset = 0;
447 tagged = NULL;
448 if (n % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"blocksize %D does not divide vector length %D", bs, n);
449 n /= bs;
450 for (i = 0; i < 2; i++) {
451 if (i) {
452 ierr = PetscMalloc1(numTagged,&tagged);CHKERRQ(ierr);
453 }
454 for (j = 0; j < n; j++) {
455
456 for (k = 0; k < numBoxes; k++) {
457 for (b = 0; b < bs; b++) {
458 PetscScalar val = vecArray[j * bs + b];
459 PetscInt l = k * bs + b;
460 VecTaggerBox box;
461 PetscBool in;
462
463 box = boxes[l];
464 #if !defined(PETSC_USE_COMPLEX)
465 in = (PetscBool) ((box.min <= val) && (val <= box.max));
466 #else
467 in = (PetscBool) ((PetscRealPart (box.min) <= PetscRealPart (val)) &&
468 (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) &&
469 (PetscRealPart (val) <= PetscRealPart (box.max)) &&
470 (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
471 #endif
472 if (!in) break;
473 }
474 if (b == bs) break;
475 }
476 if ((PetscBool)(k < numBoxes) ^ invert) {
477 if (!i) numTagged++;
478 else tagged[offset++] = j;
479 }
480 }
481 }
482 ierr = VecRestoreArrayRead(vec, &vecArray);CHKERRQ(ierr);
483 ierr = PetscFree(boxes);CHKERRQ(ierr);
484 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)vec),numTagged,tagged,PETSC_OWN_POINTER,is);CHKERRQ(ierr);
485 ierr = ISSort(*is);CHKERRQ(ierr);
486 PetscFunctionReturn(0);
487 }
488