Actual source code: swarm.c
petsc-3.8.4 2018-03-24
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petscviewer.h>
4: #include <petscdraw.h>
5: #include "data_bucket.h"
7: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
8: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
9: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
11: const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
12: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
13: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
14: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", 0 };
16: const char DMSwarmField_pid[] = "DMSwarm_pid";
17: const char DMSwarmField_rank[] = "DMSwarm_rank";
18: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
19: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
21: /*@C
22: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
23: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
25: Collective on DM
27: Input parameters:
28: + dm - a DMSwarm
29: - fieldname - the textual name given to a registered field
31: Level: beginner
33: Notes:
34: The field with name fieldname must be defined as having a data type of PetscScalar.
35: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
36: Mutiple calls to DMSwarmVectorDefineField() are permitted.
38: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
39: @*/
40: PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
41: {
42: DM_Swarm *swarm = (DM_Swarm*)dm->data;
44: PetscInt bs,n;
45: PetscScalar *array;
46: PetscDataType type;
49: if (!swarm->issetup) { DMSetUp(dm); }
50: DataBucketGetSizes(swarm->db,&n,NULL,NULL);
51: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
53: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
54: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
55: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
56: swarm->vec_field_set = PETSC_TRUE;
57: swarm->vec_field_bs = bs;
58: swarm->vec_field_nlocal = n;
59: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
60: return(0);
61: }
63: /* requires DMSwarmDefineFieldVector has been called */
64: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
65: {
66: DM_Swarm *swarm = (DM_Swarm*)dm->data;
68: Vec x;
69: char name[PETSC_MAX_PATH_LEN];
72: if (!swarm->issetup) { DMSetUp(dm); }
73: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
74: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
76: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
77: VecCreate(PetscObjectComm((PetscObject)dm),&x);
78: PetscObjectSetName((PetscObject)x,name);
79: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
80: VecSetBlockSize(x,swarm->vec_field_bs);
81: VecSetFromOptions(x);
82: *vec = x;
83: return(0);
84: }
86: /* requires DMSwarmDefineFieldVector has been called */
87: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
88: {
89: DM_Swarm *swarm = (DM_Swarm*)dm->data;
91: Vec x;
92: char name[PETSC_MAX_PATH_LEN];
95: if (!swarm->issetup) { DMSetUp(dm); }
96: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
97: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
99: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
100: VecCreate(PETSC_COMM_SELF,&x);
101: PetscObjectSetName((PetscObject)x,name);
102: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
103: VecSetBlockSize(x,swarm->vec_field_bs);
104: VecSetFromOptions(x);
105: *vec = x;
106: return(0);
107: }
109: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
110: {
111: DM_Swarm *swarm = (DM_Swarm *) dm->data;
112: DataField gfield;
113: void (*fptr)(void);
114: PetscInt bs, nlocal;
115: char name[PETSC_MAX_PATH_LEN];
119: VecGetLocalSize(*vec, &nlocal);
120: VecGetBlockSize(*vec, &bs);
121: if (nlocal/bs != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid"); /* Stale data */
122: DataBucketGetDataFieldByName(swarm->db, fieldname, &gfield);
123: /* check vector is an inplace array */
124: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
125: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
126: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
127: DataFieldRestoreAccess(gfield);
128: VecDestroy(vec);
129: return(0);
130: }
132: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
133: {
134: DM_Swarm *swarm = (DM_Swarm *) dm->data;
135: PetscDataType type;
136: PetscScalar *array;
137: PetscInt bs, n;
138: char name[PETSC_MAX_PATH_LEN];
139: PetscMPIInt commsize;
143: if (!swarm->issetup) {DMSetUp(dm);}
144: DataBucketGetSizes(swarm->db, &n, NULL, NULL);
145: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
146: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
148: MPI_Comm_size(comm, &commsize);
149: if (commsize == 1) {
150: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
151: } else {
152: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
153: }
154: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
155: PetscObjectSetName((PetscObject) *vec, name);
157: /* Set guard */
158: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
159: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
160: return(0);
161: }
163: /*@C
164: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
166: Collective on DM
168: Input parameters:
169: + dm - a DMSwarm
170: - fieldname - the textual name given to a registered field
172: Output parameter:
173: . vec - the vector
175: Level: beginner
177: Notes:
178: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
180: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
181: @*/
182: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
183: {
184: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
188: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
189: return(0);
190: }
192: /*@C
193: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
195: Collective on DM
197: Input parameters:
198: + dm - a DMSwarm
199: - fieldname - the textual name given to a registered field
201: Output parameter:
202: . vec - the vector
204: Level: beginner
206: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
207: @*/
208: PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
209: {
213: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
214: return(0);
215: }
217: /*@C
218: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
220: Collective on DM
222: Input parameters:
223: + dm - a DMSwarm
224: - fieldname - the textual name given to a registered field
226: Output parameter:
227: . vec - the vector
229: Level: beginner
231: Notes:
232: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
234: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
235: @*/
236: PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
237: {
238: MPI_Comm comm = PETSC_COMM_SELF;
242: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
243: return(0);
244: }
246: /*@C
247: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
249: Collective on DM
251: Input parameters:
252: + dm - a DMSwarm
253: - fieldname - the textual name given to a registered field
255: Output parameter:
256: . vec - the vector
258: Level: beginner
260: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
261: @*/
262: PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
263: {
267: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
268: return(0);
269: }
271: /*
272: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
273: {
274: return(0);
275: }
277: PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
278: {
279: return(0);
280: }
281: */
283: /*@C
284: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
286: Collective on DM
288: Input parameter:
289: . dm - a DMSwarm
291: Level: beginner
293: Notes:
294: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
296: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
297: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
298: @*/
299: PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
300: {
301: DM_Swarm *swarm = (DM_Swarm *) dm->data;
305: if (!swarm->field_registration_initialized) {
306: swarm->field_registration_initialized = PETSC_TRUE;
307: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_LONG); /* unique identifer */
308: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
309: }
310: return(0);
311: }
313: /*@C
314: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
316: Collective on DM
318: Input parameter:
319: . dm - a DMSwarm
321: Level: beginner
323: Notes:
324: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
326: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
327: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
328: @*/
329: PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
330: {
331: DM_Swarm *swarm = (DM_Swarm*)dm->data;
335: if (!swarm->field_registration_finalized) {
336: DataBucketFinalize(swarm->db);
337: }
338: swarm->field_registration_finalized = PETSC_TRUE;
339: return(0);
340: }
342: /*@C
343: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
345: Not collective
347: Input parameters:
348: + dm - a DMSwarm
349: . nlocal - the length of each registered field
350: - buffer - the length of the buffer used to efficient dynamic re-sizing
352: Level: beginner
354: .seealso: DMSwarmGetLocalSize()
355: @*/
356: PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
357: {
358: DM_Swarm *swarm = (DM_Swarm*)dm->data;
362: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
363: DataBucketSetSizes(swarm->db,nlocal,buffer);
364: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
365: return(0);
366: }
368: /*@C
369: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
371: Collective on DM
373: Input parameters:
374: + dm - a DMSwarm
375: - dmcell - the DM to attach to the DMSwarm
377: Level: beginner
379: Notes:
380: The attached DM (dmcell) will be queried for point location and
381: neighbor MPI-rank information if DMSwarmMigrate() is called.
383: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
384: @*/
385: PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
386: {
387: DM_Swarm *swarm = (DM_Swarm*)dm->data;
390: swarm->dmcell = dmcell;
391: return(0);
392: }
394: /*@C
395: DMSwarmGetCellDM - Fetches the attached cell DM
397: Collective on DM
399: Input parameter:
400: . dm - a DMSwarm
402: Output parameter:
403: . dmcell - the DM which was attached to the DMSwarm
405: Level: beginner
407: .seealso: DMSwarmSetCellDM()
408: @*/
409: PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
410: {
411: DM_Swarm *swarm = (DM_Swarm*)dm->data;
414: *dmcell = swarm->dmcell;
415: return(0);
416: }
418: /*@C
419: DMSwarmGetLocalSize - Retrives the local length of fields registered
421: Not collective
423: Input parameter:
424: . dm - a DMSwarm
426: Output parameter:
427: . nlocal - the length of each registered field
429: Level: beginner
431: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
432: @*/
433: PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
434: {
435: DM_Swarm *swarm = (DM_Swarm*)dm->data;
439: if (nlocal) {DataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
440: return(0);
441: }
443: /*@C
444: DMSwarmGetSize - Retrives the total length of fields registered
446: Collective on DM
448: Input parameter:
449: . dm - a DMSwarm
451: Output parameter:
452: . n - the total length of each registered field
454: Level: beginner
456: Note:
457: This calls MPI_Allreduce upon each call (inefficient but safe)
459: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
460: @*/
461: PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
462: {
463: DM_Swarm *swarm = (DM_Swarm*)dm->data;
465: PetscInt nlocal,ng;
468: DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
469: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
470: if (n) { *n = ng; }
471: return(0);
472: }
474: /*@C
475: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
477: Collective on DM
479: Input parameters:
480: + dm - a DMSwarm
481: . fieldname - the textual name to identify this field
482: . blocksize - the number of each data type
483: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
485: Level: beginner
487: Notes:
488: The textual name for each registered field must be unique.
490: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
491: @*/
492: PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
493: {
495: DM_Swarm *swarm = (DM_Swarm*)dm->data;
496: size_t size;
499: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
500: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
502: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
503: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
504: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
505: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
506: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
508: PetscDataTypeGetSize(type, &size);
509: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
510: DataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
511: {
512: DataField gfield;
514: DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);
515: DataFieldSetBlockSize(gfield,blocksize);
516: }
517: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
518: return(0);
519: }
521: /*@C
522: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
524: Collective on DM
526: Input parameters:
527: + dm - a DMSwarm
528: . fieldname - the textual name to identify this field
529: - size - the size in bytes of the user struct of each data type
531: Level: beginner
533: Notes:
534: The textual name for each registered field must be unique.
536: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
537: @*/
538: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
539: {
541: DM_Swarm *swarm = (DM_Swarm*)dm->data;
544: DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
545: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
546: return(0);
547: }
549: /*@C
550: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
552: Collective on DM
554: Input parameters:
555: + dm - a DMSwarm
556: . fieldname - the textual name to identify this field
557: . size - the size in bytes of the user data type
558: - blocksize - the number of each data type
560: Level: beginner
562: Notes:
563: The textual name for each registered field must be unique.
565: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
566: @*/
567: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
568: {
569: DM_Swarm *swarm = (DM_Swarm*)dm->data;
573: DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
574: {
575: DataField gfield;
577: DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);
578: DataFieldSetBlockSize(gfield,blocksize);
579: }
580: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
581: return(0);
582: }
584: /*@C
585: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
587: Not collective
589: Input parameters:
590: + dm - a DMSwarm
591: - fieldname - the textual name to identify this field
593: Output parameters:
594: + blocksize - the number of each data type
595: . type - the data type
596: - data - pointer to raw array
598: Level: beginner
600: Notes:
601: The array must be returned using a matching call to DMSwarmRestoreField().
603: .seealso: DMSwarmRestoreField()
604: @*/
605: PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
606: {
607: DM_Swarm *swarm = (DM_Swarm*)dm->data;
608: DataField gfield;
612: if (!swarm->issetup) { DMSetUp(dm); }
613: DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);
614: DataFieldGetAccess(gfield);
615: DataFieldGetEntries(gfield,data);
616: if (blocksize) {*blocksize = gfield->bs; }
617: if (type) { *type = gfield->petsc_type; }
618: return(0);
619: }
621: /*@C
622: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
624: Not collective
626: Input parameters:
627: + dm - a DMSwarm
628: - fieldname - the textual name to identify this field
630: Output parameters:
631: + blocksize - the number of each data type
632: . type - the data type
633: - data - pointer to raw array
635: Level: beginner
637: Notes:
638: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
640: .seealso: DMSwarmGetField()
641: @*/
642: PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
643: {
644: DM_Swarm *swarm = (DM_Swarm*)dm->data;
645: DataField gfield;
649: DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);
650: DataFieldRestoreAccess(gfield);
651: if (data) *data = NULL;
652: return(0);
653: }
655: /*@C
656: DMSwarmAddPoint - Add space for one new point in the DMSwarm
658: Not collective
660: Input parameter:
661: . dm - a DMSwarm
663: Level: beginner
665: Notes:
666: The new point will have all fields initialized to zero.
668: .seealso: DMSwarmAddNPoints()
669: @*/
670: PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
671: {
672: DM_Swarm *swarm = (DM_Swarm*)dm->data;
676: if (!swarm->issetup) {DMSetUp(dm);}
677: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
678: DataBucketAddPoint(swarm->db);
679: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
680: return(0);
681: }
683: /*@C
684: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
686: Not collective
688: Input parameters:
689: + dm - a DMSwarm
690: - npoints - the number of new points to add
692: Level: beginner
694: Notes:
695: The new point will have all fields initialized to zero.
697: .seealso: DMSwarmAddPoint()
698: @*/
699: PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
700: {
701: DM_Swarm *swarm = (DM_Swarm*)dm->data;
703: PetscInt nlocal;
706: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
707: DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
708: nlocal = nlocal + npoints;
709: DataBucketSetSizes(swarm->db,nlocal,DATA_BUCKET_BUFFER_DEFAULT);
710: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
711: return(0);
712: }
714: /*@C
715: DMSwarmRemovePoint - Remove the last point from the DMSwarm
717: Not collective
719: Input parameter:
720: . dm - a DMSwarm
722: Level: beginner
724: .seealso: DMSwarmRemovePointAtIndex()
725: @*/
726: PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
727: {
728: DM_Swarm *swarm = (DM_Swarm*)dm->data;
732: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
733: DataBucketRemovePoint(swarm->db);
734: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
735: return(0);
736: }
738: /*@C
739: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
741: Not collective
743: Input parameters:
744: + dm - a DMSwarm
745: - idx - index of point to remove
747: Level: beginner
749: .seealso: DMSwarmRemovePoint()
750: @*/
751: PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
752: {
753: DM_Swarm *swarm = (DM_Swarm*)dm->data;
757: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
758: DataBucketRemovePointAtIndex(swarm->db,idx);
759: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
760: return(0);
761: }
763: /*@C
764: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
765:
766: Not collective
767:
768: Input parameters:
769: + dm - a DMSwarm
770: . pi - the index of the point to copy
771: - pj - the point index where the copy should be located
772:
773: Level: beginner
774:
775: .seealso: DMSwarmRemovePoint()
776: @*/
777: PETSC_EXTERN PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
778: {
779: DM_Swarm *swarm = (DM_Swarm*)dm->data;
781:
783: if (!swarm->issetup) {DMSetUp(dm);}
784: DataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
785: return(0);
786: }
788: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
789: {
793: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
794: return(0);
795: }
797: /*@C
798: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
800: Collective on DM
802: Input parameters:
803: + dm - the DMSwarm
804: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
806: Notes:
807: The DM will be modified to accomodate received points.
808: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
809: Different styles of migration are supported. See DMSwarmSetMigrateType().
811: Level: advanced
813: .seealso: DMSwarmSetMigrateType()
814: @*/
815: PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
816: {
817: DM_Swarm *swarm = (DM_Swarm*)dm->data;
821: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
822: switch (swarm->migrate_type) {
823: case DMSWARM_MIGRATE_BASIC:
824: DMSwarmMigrate_Basic(dm,remove_sent_points);
825: break;
826: case DMSWARM_MIGRATE_DMCELLNSCATTER:
827: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
828: break;
829: case DMSWARM_MIGRATE_DMCELLEXACT:
830: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
831: /*DMSwarmMigrate_CellDMExact(dm,remove_sent_points);*/
832: break;
833: case DMSWARM_MIGRATE_USER:
834: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
835: /*swarm->migrate(dm,remove_sent_points);*/
836: break;
837: default:
838: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
839: break;
840: }
841: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
842: return(0);
843: }
845: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
847: /*
848: DMSwarmCollectViewCreate
850: * Applies a collection method and gathers point neighbour points into dm
852: Notes:
853: Users should call DMSwarmCollectViewDestroy() after
854: they have finished computations associated with the collected points
855: */
857: /*@C
858: DMSwarmCollectViewCreate - Applies a collection method and gathers points
859: in neighbour MPI-ranks into the DMSwarm
861: Collective on DM
863: Input parameter:
864: . dm - the DMSwarm
866: Notes:
867: Users should call DMSwarmCollectViewDestroy() after
868: they have finished computations associated with the collected points
869: Different collect methods are supported. See DMSwarmSetCollectType().
871: Level: advanced
873: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
874: @*/
875: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
876: {
878: DM_Swarm *swarm = (DM_Swarm*)dm->data;
879: PetscInt ng;
882: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
883: DMSwarmGetLocalSize(dm,&ng);
884: switch (swarm->collect_type) {
886: case DMSWARM_COLLECT_BASIC:
887: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
888: break;
889: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
890: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
891: /*DMSwarmCollect_DMDABoundingBox(dm,&ng);*/
892: break;
893: case DMSWARM_COLLECT_GENERAL:
894: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
895: /*DMSwarmCollect_General(dm,..,,..,&ng);*/
896: break;
897: default:
898: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
899: break;
900: }
901: swarm->collect_view_active = PETSC_TRUE;
902: swarm->collect_view_reset_nlocal = ng;
903: return(0);
904: }
906: /*@C
907: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
909: Collective on DM
911: Input parameters:
912: . dm - the DMSwarm
914: Notes:
915: Users should call DMSwarmCollectViewCreate() before this function is called.
917: Level: advanced
919: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
920: @*/
921: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
922: {
924: DM_Swarm *swarm = (DM_Swarm*)dm->data;
927: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
928: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
929: swarm->collect_view_active = PETSC_FALSE;
930: return(0);
931: }
933: PetscErrorCode DMSwarmSetUpPIC(DM dm)
934: {
935: PetscInt dim;
939: DMGetDimension(dm,&dim);
940: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
941: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
942: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
943: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
944: return(0);
945: }
947: /*@C
948: DMSwarmSetType - Set particular flavor of DMSwarm
950: Collective on DM
952: Input parameters:
953: + dm - the DMSwarm
954: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
956: Level: advanced
958: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
959: @*/
960: PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
961: {
962: DM_Swarm *swarm = (DM_Swarm*)dm->data;
966: swarm->swarm_type = stype;
967: if (swarm->swarm_type == DMSWARM_PIC) {
968: DMSwarmSetUpPIC(dm);
969: }
970: return(0);
971: }
973: PetscErrorCode DMSetup_Swarm(DM dm)
974: {
975: DM_Swarm *swarm = (DM_Swarm*)dm->data;
977: PetscMPIInt rank;
978: PetscInt p,npoints,*rankval;
981: if (swarm->issetup) return(0);
983: swarm->issetup = PETSC_TRUE;
985: if (swarm->swarm_type == DMSWARM_PIC) {
986: /* check dmcell exists */
987: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
989: if (swarm->dmcell->ops->locatepointssubdomain) {
990: /* check methods exists for exact ownership identificiation */
991: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
992: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
993: } else {
994: /* check methods exist for point location AND rank neighbor identification */
995: if (swarm->dmcell->ops->locatepoints) {
996: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->LocatePoints\n");
997: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
999: if (swarm->dmcell->ops->getneighbors) {
1000: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1001: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1003: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1004: }
1005: }
1007: DMSwarmFinalizeFieldRegister(dm);
1009: /* check some fields were registered */
1010: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1012: /* check local sizes were set */
1013: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1015: /* initialize values in pid and rank placeholders */
1016: /* TODO: [pid - use MPI_Scan] */
1017: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1018: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1019: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1020: for (p=0; p<npoints; p++) {
1021: rankval[p] = (PetscInt)rank;
1022: }
1023: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1024: return(0);
1025: }
1027: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1029: PetscErrorCode DMDestroy_Swarm(DM dm)
1030: {
1031: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1035: DataBucketDestroy(&swarm->db);
1036: if (swarm->sort_context) {
1037: DMSwarmSortDestroy(&swarm->sort_context);
1038: }
1039: PetscFree(swarm);
1040: return(0);
1041: }
1043: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1044: {
1045: DM cdm;
1046: PetscDraw draw;
1047: PetscReal *coords, oldPause;
1048: PetscInt Np, p, bs;
1052: PetscViewerDrawGetDraw(viewer, 0, &draw);
1053: DMSwarmGetCellDM(dm, &cdm);
1054: PetscDrawGetPause(draw, &oldPause);
1055: PetscDrawSetPause(draw, 0.0);
1056: DMView(cdm, viewer);
1057: PetscDrawSetPause(draw, oldPause);
1059: DMSwarmGetLocalSize(dm, &Np);
1060: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1061: for (p = 0; p < Np; ++p) {
1062: const PetscInt i = p*bs;
1064: PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);
1065: }
1066: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1067: PetscDrawFlush(draw);
1068: PetscDrawPause(draw);
1069: PetscDrawSave(draw);
1070: return(0);
1071: }
1073: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1074: {
1075: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1076: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1082: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1083: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1084: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1085: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1086: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1087: if (iascii) {
1088: DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1089: } else if (ibinary) {
1090: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1091: } else if (ishdf5) {
1092: #if defined(PETSC_HAVE_HDF5)
1093: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1094: #else
1095: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1096: #endif
1097: } else if (isvtk) {
1098: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1099: } else if (isdraw) {
1100: DMSwarmView_Draw(dm, viewer);
1101: }
1102: return(0);
1103: }
1105: /*MC
1107: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1108: This implementation was designed for particle methods in which the underlying
1109: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1111: User data can be represented by DMSwarm through a registering "fields".
1112: To register a field, the user must provide;
1113: (a) a unique name;
1114: (b) the data type (or size in bytes);
1115: (c) the block size of the data.
1117: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1118: on a set of of particles. Then the following code could be used
1120: $ DMSwarmInitializeFieldRegister(dm)
1121: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1122: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1123: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1124: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1125: $ DMSwarmFinalizeFieldRegister(dm)
1127: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1128: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1130: To support particle methods, "migration" techniques are provided. These methods migrate data
1131: between MPI-ranks.
1133: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1134: As a DMSwarm may internally define and store values of different data types,
1135: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1136: fields should be used to define a Vec object via
1137: DMSwarmVectorDefineField()
1138: The specified field can can changed be changed at any time - thereby permitting vectors
1139: compatable with different fields to be created.
1141: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1142: DMSwarmCreateGlobalVectorFromField()
1143: Here the data defining the field in the DMSwarm is shared with a Vec.
1144: This is inherently unsafe if you alter the size of the field at any time between
1145: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1146: If the local size of the DMSwarm does not match the local size of the global vector
1147: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1149: Additional high-level support is provided for Particle-In-Cell methods.
1150: Please refer to the man page for DMSwarmSetType().
1151:
1152: Level: beginner
1154: .seealso: DMType, DMCreate(), DMSetType()
1155: M*/
1156: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1157: {
1158: DM_Swarm *swarm;
1163: PetscNewLog(dm,&swarm);
1164: dm->data = swarm;
1166: DataBucketCreate(&swarm->db);
1167: DMSwarmInitializeFieldRegister(dm);
1169: swarm->vec_field_set = PETSC_FALSE;
1170: swarm->issetup = PETSC_FALSE;
1171: swarm->swarm_type = DMSWARM_BASIC;
1172: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1173: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1174: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1176: swarm->dmcell = NULL;
1177: swarm->collect_view_active = PETSC_FALSE;
1178: swarm->collect_view_reset_nlocal = -1;
1180: dm->dim = 0;
1181: dm->ops->view = DMView_Swarm;
1182: dm->ops->load = NULL;
1183: dm->ops->setfromoptions = NULL;
1184: dm->ops->clone = NULL;
1185: dm->ops->setup = DMSetup_Swarm;
1186: dm->ops->createdefaultsection = NULL;
1187: dm->ops->createdefaultconstraints = NULL;
1188: dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1189: dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1190: dm->ops->getlocaltoglobalmapping = NULL;
1191: dm->ops->createfieldis = NULL;
1192: dm->ops->createcoordinatedm = NULL;
1193: dm->ops->getcoloring = NULL;
1194: dm->ops->creatematrix = NULL;
1195: dm->ops->createinterpolation = NULL;
1196: dm->ops->getaggregates = NULL;
1197: dm->ops->getinjection = NULL;
1198: dm->ops->refine = NULL;
1199: dm->ops->coarsen = NULL;
1200: dm->ops->refinehierarchy = NULL;
1201: dm->ops->coarsenhierarchy = NULL;
1202: dm->ops->globaltolocalbegin = NULL;
1203: dm->ops->globaltolocalend = NULL;
1204: dm->ops->localtoglobalbegin = NULL;
1205: dm->ops->localtoglobalend = NULL;
1206: dm->ops->destroy = DMDestroy_Swarm;
1207: dm->ops->createsubdm = NULL;
1208: dm->ops->getdimpoints = NULL;
1209: dm->ops->locatepoints = NULL;
1210: return(0);
1211: }