1 2 #include "data_bucket.h" 3 4 /* string helpers */ 5 PetscErrorCode StringInList( const char name[], const PetscInt N, const DataField gfield[], PetscBool *val ) 6 { 7 PetscInt i; 8 9 *val = PETSC_FALSE; 10 for( i=0; i<N; i++ ) { 11 if( strcmp( name, gfield[i]->name ) == 0 ) { 12 *val = PETSC_TRUE; 13 PetscFunctionReturn(0); 14 } 15 } 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode StringFindInList( const char name[], const PetscInt N, const DataField gfield[], PetscInt *index ) 20 { 21 PetscInt i; 22 23 *index = -1; 24 for( i=0; i<N; i++ ) { 25 if( strcmp( name, gfield[i]->name ) == 0 ) { 26 *index = i; 27 PetscFunctionReturn(0); 28 } 29 } 30 PetscFunctionReturn(0); 31 } 32 33 #undef __FUNCT__ 34 #define __FUNCT__ "DataFieldCreate" 35 PetscErrorCode DataFieldCreate( const char registeration_function[], const char name[], const size_t size, const PetscInt L, DataField *DF ) 36 { 37 DataField df; 38 39 df = malloc( sizeof(struct _p_DataField) ); 40 memset( df, 0, sizeof(struct _p_DataField) ); 41 42 43 asprintf( &df->registeration_function, "%s", registeration_function ); 44 asprintf( &df->name, "%s", name ); 45 df->atomic_size = size; 46 df->L = L; 47 48 df->data = malloc( size * L ); /* allocate something so we don't have to reallocate */ 49 memset( df->data, 0, size * L ); 50 51 *DF = df; 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "DataFieldDestroy" 57 PetscErrorCode DataFieldDestroy( DataField *DF ) 58 { 59 DataField df = *DF; 60 61 free( df->registeration_function ); 62 free( df->name ); 63 free( df->data ); 64 free(df); 65 66 *DF = NULL; 67 PetscFunctionReturn(0); 68 } 69 70 /* data bucket */ 71 #undef __FUNCT__ 72 #define __FUNCT__ "DataBucketCreate" 73 PetscErrorCode DataBucketCreate( DataBucket *DB ) 74 { 75 DataBucket db; 76 77 78 db = malloc( sizeof(struct _p_DataBucket) ); 79 memset( db, 0, sizeof(struct _p_DataBucket) ); 80 81 db->finalised = PETSC_FALSE; 82 83 /* create empty spaces for fields */ 84 db->L = 0; 85 db->buffer = 1; 86 db->allocated = 1; 87 88 db->nfields = 0; 89 db->field = malloc(sizeof(DataField)); 90 91 *DB = db; 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "DataBucketDestroy" 97 PetscErrorCode DataBucketDestroy( DataBucket *DB ) 98 { 99 DataBucket db = *DB; 100 PetscInt f; 101 102 /* release fields */ 103 for( f=0; f<db->nfields; f++ ) { 104 DataFieldDestroy(&db->field[f]); 105 } 106 107 /* this will catch the initially allocated objects in the event that no fields are registered */ 108 if(db->field!=NULL) { 109 free(db->field); 110 } 111 112 free(db); 113 114 *DB = NULL; 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "DataBucketRegisterField" 120 PetscErrorCode DataBucketRegisterField( 121 DataBucket db, 122 const char registeration_function[], 123 const char field_name[], 124 size_t atomic_size, DataField *_gfield ) 125 { 126 PetscBool val; 127 DataField *field,fp; 128 129 /* check we haven't finalised the registration of fields */ 130 /* 131 if(db->finalised==PETSC_TRUE) { 132 printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n"); 133 ERROR(); 134 } 135 */ 136 137 /* check for repeated name */ 138 StringInList( field_name, db->nfields, (const DataField*)db->field, &val ); 139 if(val == PETSC_TRUE ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field %s already exists. Cannot add same field twice",field_name); 140 141 /* create new space for data */ 142 field = realloc( db->field, sizeof(DataField)*(db->nfields+1)); 143 db->field = field; 144 145 /* add field */ 146 DataFieldCreate( registeration_function, field_name, atomic_size, db->allocated, &fp ); 147 db->field[ db->nfields ] = fp; 148 149 db->nfields++; 150 151 if(_gfield!=NULL){ 152 *_gfield = fp; 153 } 154 PetscFunctionReturn(0); 155 } 156 157 /* 158 #define DataBucketRegisterField(db,name,size,k) {\ 159 char *location;\ 160 asprintf(&location,"Registered by %s() at line %d within file %s", __FUNCTION__, __LINE__, __FILE__);\ 161 _DataBucketRegisterField( (db), location, (name), (size), (k) );\ 162 free(location);\ 163 } 164 */ 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "DataBucketGetDataFieldByName" 168 PetscErrorCode DataBucketGetDataFieldByName(DataBucket db,const char name[],DataField *gfield) 169 { 170 PetscInt idx; 171 PetscBool found; 172 173 StringInList(name,db->nfields,(const DataField*)db->field,&found); 174 if (!found) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot find DataField with name %s",name); 175 176 StringFindInList(name,db->nfields,(const DataField*)db->field,&idx); 177 178 *gfield = db->field[idx]; 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "DataBucketQueryDataFieldByName" 184 PetscErrorCode DataBucketQueryDataFieldByName(DataBucket db,const char name[],PetscBool *found) 185 { 186 *found = PETSC_FALSE; 187 StringInList(name,db->nfields,(const DataField*)db->field,found); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "DataBucketFinalize" 193 PetscErrorCode DataBucketFinalize(DataBucket db) 194 { 195 db->finalised = PETSC_TRUE; 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "DataFieldGetNumEntries" 201 PetscErrorCode DataFieldGetNumEntries(DataField df, PetscInt *sum) 202 { 203 *sum = df->L; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DataFieldSetSize" 209 PetscErrorCode DataFieldSetSize( DataField df, const PetscInt new_L ) 210 { 211 void *tmp_data; 212 213 if (new_L <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot set size of DataField to be <= 0"); 214 215 if (new_L == df->L) PetscFunctionReturn(0); 216 217 if (new_L > df->L) { 218 219 tmp_data = realloc( df->data, df->atomic_size * (new_L) ); 220 df->data = tmp_data; 221 222 /* init new contents */ 223 memset( ( ((char*)df->data)+df->L*df->atomic_size), 0, (new_L-df->L)*df->atomic_size ); 224 225 } else { 226 /* reallocate pointer list, add +1 in case new_L = 0 */ 227 tmp_data = realloc( df->data, df->atomic_size * (new_L+1) ); 228 df->data = tmp_data; 229 } 230 231 df->L = new_L; 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "DataFieldZeroBlock" 237 PetscErrorCode DataFieldZeroBlock( DataField df, const PetscInt start, const PetscInt end ) 238 { 239 if (start > end) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) > end(%D)",start,end); 240 241 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) < 0",start); 242 243 if (end > df->L) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if end(%D) >= array size(%D)",end,df->L); 244 245 memset( ( ((char*)df->data)+start*df->atomic_size), 0, (end-start)*df->atomic_size ); 246 PetscFunctionReturn(0); 247 } 248 249 /* 250 A negative buffer value will simply be ignored and the old buffer value will be used. 251 */ 252 #undef __FUNCT__ 253 #define __FUNCT__ "DataBucketSetSizes" 254 PetscErrorCode DataBucketSetSizes( DataBucket db, const PetscInt L, const PetscInt buffer ) 255 { 256 PetscInt current_allocated,new_used,new_unused,new_buffer,new_allocated,f; 257 258 if (db->finalised == PETSC_FALSE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"You must call DataBucketFinalize() before DataBucketSetSizes()"); 259 260 current_allocated = db->allocated; 261 262 new_used = L; 263 new_unused = current_allocated - new_used; 264 new_buffer = db->buffer; 265 if( buffer >= 0 ) { /* update the buffer value */ 266 new_buffer = buffer; 267 } 268 new_allocated = new_used + new_buffer; 269 270 /* action */ 271 if ( new_allocated > current_allocated ) { 272 /* increase size to new_used + new_buffer */ 273 for( f=0; f<db->nfields; f++ ) { 274 DataFieldSetSize( db->field[f], new_allocated ); 275 } 276 277 db->L = new_used; 278 db->buffer = new_buffer; 279 db->allocated = new_used + new_buffer; 280 } 281 else { 282 if( new_unused > 2 * new_buffer ) { 283 284 /* shrink array to new_used + new_buffer */ 285 for( f=0; f<db->nfields; f++ ) { 286 DataFieldSetSize( db->field[f], new_allocated ); 287 } 288 289 db->L = new_used; 290 db->buffer = new_buffer; 291 db->allocated = new_used + new_buffer; 292 } 293 else { 294 db->L = new_used; 295 db->buffer = new_buffer; 296 } 297 } 298 299 /* zero all entries from db->L to db->allocated */ 300 for( f=0; f<db->nfields; f++ ) { 301 DataField field = db->field[f]; 302 DataFieldZeroBlock(field, db->L,db->allocated); 303 } 304 PetscFunctionReturn(0); 305 } 306 307 #undef __FUNCT__ 308 #define __FUNCT__ "DataBucketSetInitialSizes" 309 PetscErrorCode DataBucketSetInitialSizes( DataBucket db, const PetscInt L, const PetscInt buffer ) 310 { 311 PetscInt f; 312 DataBucketSetSizes(db,L,buffer); 313 314 for( f=0; f<db->nfields; f++ ) { 315 DataField field = db->field[f]; 316 DataFieldZeroBlock(field,0,db->allocated); 317 } 318 PetscFunctionReturn(0); 319 } 320 321 #undef __FUNCT__ 322 #define __FUNCT__ "DataBucketGetSizes" 323 PetscErrorCode DataBucketGetSizes( DataBucket db, PetscInt *L, PetscInt *buffer, PetscInt *allocated ) 324 { 325 if (L) { *L = db->L; } 326 if (buffer) { *buffer = db->buffer; } 327 if (allocated) { *allocated = db->allocated; } 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "DataBucketGetGlobalSizes" 333 PetscErrorCode DataBucketGetGlobalSizes(MPI_Comm comm, DataBucket db, PetscInt *L, PetscInt *buffer, PetscInt *allocated ) 334 { 335 PetscInt _L,_buffer,_allocated; 336 PetscInt ierr; 337 338 _L = db->L; 339 _buffer = db->buffer; 340 _allocated = db->allocated; 341 342 if (L) { ierr = MPI_Allreduce(&_L,L,1,MPIU_INT,MPI_SUM,comm); } 343 if (buffer) { ierr = MPI_Allreduce(&_buffer,buffer,1,MPIU_INT,MPI_SUM,comm); } 344 if (allocated) { ierr = MPI_Allreduce(&_allocated,allocated,1,MPIU_INT,MPI_SUM,comm); } 345 PetscFunctionReturn(0); 346 } 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "DataBucketGetGlobalSizes" 350 PetscErrorCode DataBucketGetDataFields( DataBucket db, PetscInt *L, DataField *fields[] ) 351 { 352 if (L) { *L = db->nfields; } 353 if (fields) { *fields = db->field; } 354 PetscFunctionReturn(0); 355 } 356 357 #undef __FUNCT__ 358 #define __FUNCT__ "DataFieldGetAccess" 359 PetscErrorCode DataFieldGetAccess( const DataField gfield ) 360 { 361 if (gfield->active) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is already active. You must call DataFieldRestoreAccess()",gfield->name); 362 363 gfield->active = PETSC_TRUE; 364 PetscFunctionReturn(0); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "DataFieldAccessPoint" 369 PetscErrorCode DataFieldAccessPoint( const DataField gfield, const PetscInt pid, void **ctx_p ) 370 { 371 #ifdef DATAFIELD_POINT_ACCESS_GUARD 372 /* debug mode */ 373 /* check poPetscInt is valid */ 374 if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 375 if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L); 376 377 if (gfield->active == PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess() before poPetscInt data can be retrivied",gfield->name); 378 #endif 379 380 //*ctx_p = (void*)( ((char*)gfield->data) + pid * gfield->atomic_size); 381 *ctx_p = __DATATFIELD_point_access(gfield->data,pid,gfield->atomic_size); 382 PetscFunctionReturn(0); 383 } 384 385 #undef __FUNCT__ 386 #define __FUNCT__ "DataFieldAccessPointOffset" 387 PetscErrorCode DataFieldAccessPointOffset( const DataField gfield, const size_t offset, const PetscInt pid, void **ctx_p ) 388 { 389 #ifdef DATAFIELD_POINT_ACCESS_GUARD 390 /* debug mode */ 391 392 /* check poPetscInt is valid */ 393 /* if( offset < 0 ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"offset must be >= 0");*//* Note compiler realizes this can never happen with an unsigned PetscInt */ 394 if (offset >= gfield->atomic_size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"offset must be < %zu",gfield->atomic_size); 395 396 /* check poPetscInt is valid */ 397 if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 398 if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L); 399 400 if (gfield->active==PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess() before point data can be retrivied",gfield->name); 401 #endif 402 403 *ctx_p = __DATATFIELD_point_access_offset(gfield->data,pid,gfield->atomic_size,offset); 404 PetscFunctionReturn(0); 405 } 406 407 #undef __FUNCT__ 408 #define __FUNCT__ "DataFieldAccessPointOffset" 409 PetscErrorCode DataFieldRestoreAccess( DataField gfield ) 410 { 411 if (gfield->active == PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess()", gfield->name ); 412 413 gfield->active = PETSC_FALSE; 414 PetscFunctionReturn(0); 415 } 416 417 #undef __FUNCT__ 418 #define __FUNCT__ "DataFieldVerifyAccess" 419 PetscErrorCode DataFieldVerifyAccess( const DataField gfield, const size_t size) 420 { 421 #ifdef DATAFIELD_POINT_ACCESS_GUARD 422 if (gfield->atomic_size != size) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" must be mapped to %zu bytes, your intended structure is %zu bytes in length.", 423 gfield->name, gfield->atomic_size, size ); 424 #endif 425 PetscFunctionReturn(0); 426 } 427 428 #undef __FUNCT__ 429 #define __FUNCT__ "DataFieldGetAtomicSize" 430 PetscErrorCode DataFieldGetAtomicSize(const DataField gfield,size_t *size) 431 { 432 if (size) { *size = gfield->atomic_size; } 433 PetscFunctionReturn(0); 434 } 435 436 #undef __FUNCT__ 437 #define __FUNCT__ "DataFieldGetEntries" 438 PetscErrorCode DataFieldGetEntries(const DataField gfield,void **data) 439 { 440 if (data) { 441 *data = gfield->data; 442 } 443 PetscFunctionReturn(0); 444 } 445 446 #undef __FUNCT__ 447 #define __FUNCT__ "DataFieldRestoreEntries" 448 PetscErrorCode DataFieldRestoreEntries(const DataField gfield,void **data) 449 { 450 if (data) { 451 *data = NULL; 452 } 453 PetscFunctionReturn(0); 454 } 455 456 /* y = x */ 457 #undef __FUNCT__ 458 #define __FUNCT__ "DataBucketCopyPoint" 459 PetscErrorCode DataBucketCopyPoint( const DataBucket xb, const PetscInt pid_x, 460 const DataBucket yb, const PetscInt pid_y ) 461 { 462 PetscInt f; 463 for( f=0; f<xb->nfields; f++ ) { 464 void *dest; 465 void *src; 466 467 DataFieldGetAccess( xb->field[f] ); 468 if (xb!=yb) { DataFieldGetAccess( yb->field[f] ); } 469 470 DataFieldAccessPoint( xb->field[f],pid_x, &src ); 471 DataFieldAccessPoint( yb->field[f],pid_y, &dest ); 472 473 memcpy( dest, src, xb->field[f]->atomic_size ); 474 475 DataFieldRestoreAccess( xb->field[f] ); 476 if (xb!=yb) { DataFieldRestoreAccess( yb->field[f] ); } 477 } 478 PetscFunctionReturn(0); 479 } 480 481 #undef __FUNCT__ 482 #define __FUNCT__ "DataBucketCreateFromSubset" 483 PetscErrorCode DataBucketCreateFromSubset( DataBucket DBIn, const PetscInt N, const PetscInt list[], DataBucket *DB ) 484 { 485 PetscInt nfields; 486 DataField *fields; 487 DataBucketCreate(DB); 488 PetscInt f,L,buffer,allocated,p; 489 490 /* copy contents of DBIn */ 491 DataBucketGetDataFields(DBIn,&nfields,&fields); 492 DataBucketGetSizes(DBIn,&L,&buffer,&allocated); 493 494 for(f=0;f<nfields;f++) { 495 DataBucketRegisterField(*DB,"DataBucketCreateFromSubset",fields[f]->name,fields[f]->atomic_size,NULL); 496 } 497 DataBucketFinalize(*DB); 498 499 DataBucketSetSizes(*DB,L,buffer); 500 501 /* now copy the desired guys from DBIn => DB */ 502 for( p=0; p<N; p++ ) { 503 DataBucketCopyPoint(DBIn,list[p], *DB,p); 504 } 505 PetscFunctionReturn(0); 506 } 507 508 // insert into an exisitng location 509 #undef __FUNCT__ 510 #define __FUNCT__ "DataFieldInsertPoint" 511 PetscErrorCode DataFieldInsertPoint( const DataField field, const PetscInt index, const void *ctx ) 512 { 513 514 #ifdef DATAFIELD_POINT_ACCESS_GUARD 515 /* check poPetscInt is valid */ 516 if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 517 if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L); 518 #endif 519 520 // memcpy( (void*)((char*)field->data + index*field->atomic_size), ctx, field->atomic_size ); 521 memcpy( __DATATFIELD_point_access(field->data,index,field->atomic_size), ctx, field->atomic_size ); 522 PetscFunctionReturn(0); 523 } 524 525 // remove data at index - replace with last point 526 #undef __FUNCT__ 527 #define __FUNCT__ "DataBucketRemovePointAtIndex" 528 PetscErrorCode DataBucketRemovePointAtIndex( const DataBucket db, const PetscInt index ) 529 { 530 PetscInt f; 531 532 #ifdef DATAFIELD_POINT_ACCESS_GUARD 533 /* check poPetscInt is valid */ 534 if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 535 if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->L+db->buffer); 536 #endif 537 538 if (index >= db->L) { /* this point is not in the list - no need to error, but I will anyway */ 539 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"You should not be trying to remove point at index=%D since it's < db->L = %D", index, db->L ); 540 } 541 542 if (index != db->L-1) { /* not last point in list */ 543 for( f=0; f<db->nfields; f++ ) { 544 DataField field = db->field[f]; 545 546 /* copy then remove */ 547 DataFieldCopyPoint( db->L-1,field, index,field ); 548 549 //DataFieldZeroPoint(field,index); 550 } 551 } 552 553 /* decrement size */ 554 /* this will zero out an crap at the end of the list */ 555 DataBucketRemovePoint(db); 556 PetscFunctionReturn(0); 557 } 558 559 /* copy x into y */ 560 #undef __FUNCT__ 561 #define __FUNCT__ "DataFieldCopyPoint" 562 PetscErrorCode DataFieldCopyPoint( const PetscInt pid_x, const DataField field_x, 563 const PetscInt pid_y, const DataField field_y ) 564 { 565 566 #ifdef DATAFIELD_POINT_ACCESS_GUARD 567 /* check poPetscInt is valid */ 568 if (pid_x < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be >= 0"); 569 if (pid_x >= field_x->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be < %D",field_x->L); 570 571 if (pid_y < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be >= 0"); 572 if (pid_y >= field_y->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be < %D",field_y->L); 573 574 if( field_y->atomic_size != field_x->atomic_size ) { 575 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"atomic size must match"); 576 } 577 #endif 578 /* 579 memcpy( (void*)((char*)field_y->data + pid_y*field_y->atomic_size), 580 (void*)((char*)field_x->data + pid_x*field_x->atomic_size), 581 field_x->atomic_size ); 582 */ 583 memcpy( __DATATFIELD_point_access(field_y->data,pid_y,field_y->atomic_size), 584 __DATATFIELD_point_access(field_x->data,pid_x,field_x->atomic_size), 585 field_y->atomic_size ); 586 PetscFunctionReturn(0); 587 } 588 589 590 // zero only the datafield at this point 591 #undef __FUNCT__ 592 #define __FUNCT__ "DataFieldZeroPoint" 593 PetscErrorCode DataFieldZeroPoint( const DataField field, const PetscInt index ) 594 { 595 #ifdef DATAFIELD_POINT_ACCESS_GUARD 596 /* check poPetscInt is valid */ 597 if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 598 if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L); 599 #endif 600 601 // memset( (void*)((char*)field->data + index*field->atomic_size), 0, field->atomic_size ); 602 memset( __DATATFIELD_point_access(field->data,index,field->atomic_size), 0, field->atomic_size ); 603 PetscFunctionReturn(0); 604 } 605 606 // zero ALL data for this point 607 #undef __FUNCT__ 608 #define __FUNCT__ "DataBucketZeroPoint" 609 PetscErrorCode DataBucketZeroPoint( const DataBucket db, const PetscInt index ) 610 { 611 PetscInt f; 612 613 /* check poPetscInt is valid */ 614 if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0"); 615 if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->allocated); 616 617 for(f=0;f<db->nfields;f++){ 618 DataField field = db->field[f]; 619 620 DataFieldZeroPoint(field,index); 621 } 622 PetscFunctionReturn(0); 623 } 624 625 /* increment */ 626 #undef __FUNCT__ 627 #define __FUNCT__ "DataBucketAddPoint" 628 PetscErrorCode DataBucketAddPoint( DataBucket db ) 629 { 630 DataBucketSetSizes( db, db->L+1, -1 ); 631 PetscFunctionReturn(0); 632 } 633 634 /* decrement */ 635 #undef __FUNCT__ 636 #define __FUNCT__ "DataBucketRemovePoint" 637 PetscErrorCode DataBucketRemovePoint( DataBucket db ) 638 { 639 DataBucketSetSizes( db, db->L-1, -1 ); 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "_DataFieldViewBinary" 645 PetscErrorCode _DataFieldViewBinary(DataField field, FILE *fp ) 646 { 647 fprintf(fp,"<DataField>\n"); 648 fprintf(fp,"%d\n", field->L); 649 fprintf(fp,"%zu\n",field->atomic_size); 650 fprintf(fp,"%s\n", field->registeration_function); 651 fprintf(fp,"%s\n", field->name); 652 653 fwrite(field->data, field->atomic_size, field->L, fp); 654 /* 655 printf(" ** wrote %zu bytes for DataField \"%s\" \n", field->atomic_size * field->L, field->name ); 656 */ 657 fprintf(fp,"\n</DataField>\n"); 658 PetscFunctionReturn(0); 659 } 660 661 #undef __FUNCT__ 662 #define __FUNCT__ "_DataBucketRegisterFieldFromFile" 663 PetscErrorCode _DataBucketRegisterFieldFromFile( FILE *fp, DataBucket db ) 664 { 665 PetscBool val; 666 DataField *field; 667 668 DataField gfield; 669 char dummy[100]; 670 char registeration_function[5000]; 671 char field_name[5000]; 672 PetscInt L; 673 size_t atomic_size,strL; 674 675 676 /* check we haven't finalised the registration of fields */ 677 /* 678 if(db->finalised==PETSC_TRUE) { 679 printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n"); 680 ERROR(); 681 } 682 */ 683 684 685 /* read file contents */ 686 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 687 688 fscanf( fp, "%d\n",&L); //printf("read(L): %d\n", L); 689 690 fscanf( fp, "%zu\n",&atomic_size); //printf("read(size): %zu\n",atomic_size); 691 692 fgets(registeration_function,4999,fp); //printf("read(reg func): %s", registeration_function ); 693 strL = strlen(registeration_function); 694 if(strL>1){ 695 registeration_function[strL-1] = 0; 696 } 697 698 fgets(field_name,4999,fp); //printf("read(name): %s", field_name ); 699 strL = strlen(field_name); 700 if(strL>1){ 701 field_name[strL-1] = 0; 702 } 703 704 #ifdef DATA_BUCKET_LOG 705 printf(" ** read L=%d; atomic_size=%zu; reg_func=\"%s\"; name=\"%s\" \n", L,atomic_size,registeration_function,field_name); 706 #endif 707 708 709 /* check for repeated name */ 710 StringInList( field_name, db->nfields, (const DataField*)db->field, &val ); 711 if (val == PETSC_TRUE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot add same field twice"); 712 713 /* create new space for data */ 714 field = realloc( db->field, sizeof(DataField)*(db->nfields+1)); 715 db->field = field; 716 717 /* add field */ 718 DataFieldCreate( registeration_function, field_name, atomic_size, L, &gfield ); 719 720 /* copy contents of file */ 721 fread(gfield->data, gfield->atomic_size, gfield->L, fp); 722 #ifdef DATA_BUCKET_LOG 723 printf(" ** read %zu bytes for DataField \"%s\" \n", gfield->atomic_size * gfield->L, field_name ); 724 #endif 725 /* finish reading meta data */ 726 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 727 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 728 729 db->field[ db->nfields ] = gfield; 730 731 db->nfields++; 732 PetscFunctionReturn(0); 733 } 734 735 #undef __FUNCT__ 736 #define __FUNCT__ "_DataBucketViewAscii_HeaderWrite_v00" 737 PetscErrorCode _DataBucketViewAscii_HeaderWrite_v00(FILE *fp) 738 { 739 fprintf(fp,"<DataBucketHeader>\n"); 740 fprintf(fp,"type=DataBucket\n"); 741 fprintf(fp,"format=ascii\n"); 742 fprintf(fp,"version=0.0\n"); 743 fprintf(fp,"options=\n"); 744 fprintf(fp,"</DataBucketHeader>\n"); 745 PetscFunctionReturn(0); 746 } 747 748 #undef __FUNCT__ 749 #define __FUNCT__ "_DataBucketViewAscii_HeaderRead_v00" 750 PetscErrorCode _DataBucketViewAscii_HeaderRead_v00(FILE *fp) 751 { 752 char dummy[100]; 753 size_t strL; 754 755 // header open 756 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 757 758 // type 759 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 760 strL = strlen(dummy); 761 if(strL>1) { dummy[strL-1] = 0; } 762 if(strcmp(dummy,"type=DataBucket")!=0) { 763 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Data file doesn't contain a DataBucket type"); 764 } 765 766 // format 767 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 768 769 // version 770 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 771 strL = strlen(dummy); 772 if(strL>1) { dummy[strL-1] = 0; } 773 if (strcmp(dummy,"version=0.0") != 0) { 774 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"DataBucket file must be parsed with version=0.0 : You tried %s", dummy); 775 } 776 777 // options 778 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 779 // header close 780 fgets(dummy,99,fp); //printf("read(header): %s", dummy ); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "_DataBucketLoadFromFileBinary_SEQ" 786 PetscErrorCode _DataBucketLoadFromFileBinary_SEQ(const char filename[], DataBucket *_db) 787 { 788 DataBucket db; 789 FILE *fp; 790 PetscInt L,buffer,f,nfields; 791 792 793 #ifdef DATA_BUCKET_LOG 794 printf("** DataBucketLoadFromFile **\n"); 795 #endif 796 797 /* open file */ 798 fp = fopen(filename,"rb"); 799 if (fp==NULL){ 800 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file with name %s", filename); 801 } 802 803 /* read header */ 804 _DataBucketViewAscii_HeaderRead_v00(fp); 805 806 fscanf(fp,"%d\n%d\n%d\n",&L,&buffer,&nfields); 807 808 DataBucketCreate(&db); 809 810 for( f=0; f<nfields; f++ ) { 811 _DataBucketRegisterFieldFromFile(fp,db); 812 } 813 fclose(fp); 814 815 DataBucketFinalize(db); 816 817 818 /* 819 DataBucketSetSizes(db,L,buffer); 820 */ 821 db->L = L; 822 db->buffer = buffer; 823 db->allocated = L + buffer; 824 825 *_db = db; 826 PetscFunctionReturn(0); 827 } 828 829 #undef __FUNCT__ 830 #define __FUNCT__ "DataBucketLoadFromFile" 831 PetscErrorCode DataBucketLoadFromFile(MPI_Comm comm,const char filename[], DataBucketViewType type, DataBucket *db) 832 { 833 PetscMPIInt nproc,rank; 834 835 MPI_Comm_size(comm,&nproc); 836 MPI_Comm_rank(comm,&rank); 837 838 #ifdef DATA_BUCKET_LOG 839 printf("** DataBucketLoadFromFile **\n"); 840 #endif 841 if (type == DATABUCKET_VIEW_STDOUT) { 842 843 } else if (type == DATABUCKET_VIEW_ASCII) { 844 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure"); 845 } else if (type == DATABUCKET_VIEW_BINARY) { 846 if (nproc == 1) { 847 _DataBucketLoadFromFileBinary_SEQ(filename,db); 848 } else { 849 char *name; 850 851 asprintf(&name,"%s_p%1.5d",filename, rank ); 852 _DataBucketLoadFromFileBinary_SEQ(name,db); 853 free(name); 854 } 855 } else { 856 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer requested"); 857 } 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNCT__ 862 #define __FUNCT__ "_DataBucketViewBinary" 863 PetscErrorCode _DataBucketViewBinary(DataBucket db,const char filename[]) 864 { 865 FILE *fp = NULL; 866 PetscInt f; 867 868 fp = fopen(filename,"wb"); 869 if (fp == NULL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file with name %s", filename); 870 871 /* db header */ 872 _DataBucketViewAscii_HeaderWrite_v00(fp); 873 874 /* meta-data */ 875 fprintf(fp,"%d\n%d\n%d\n", db->L,db->buffer,db->nfields); 876 877 for( f=0; f<db->nfields; f++ ) { 878 /* load datafields */ 879 _DataFieldViewBinary(db->field[f],fp); 880 } 881 882 fclose(fp); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "DataBucketView_SEQ" 888 PetscErrorCode DataBucketView_SEQ(DataBucket db,const char filename[],DataBucketViewType type) 889 { 890 switch (type) { 891 case DATABUCKET_VIEW_STDOUT: 892 { 893 PetscInt f; 894 double memory_usage_total = 0.0; 895 896 printf("DataBucketView(SEQ): (\"%s\")\n",filename); 897 printf(" L = %d \n", db->L ); 898 printf(" buffer = %d \n", db->buffer ); 899 printf(" allocated = %d \n", db->allocated ); 900 901 printf(" nfields registered = %d \n", db->nfields ); 902 for( f=0; f<db->nfields; f++ ) { 903 double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6; 904 905 printf(" [%3d]: field name ==>> %30s : Mem. usage = %1.2e (MB) \n", f, db->field[f]->name, memory_usage_f ); 906 memory_usage_total += memory_usage_f; 907 } 908 printf(" Total mem. usage = %1.2e (MB) \n", memory_usage_total ); 909 } 910 break; 911 912 case DATABUCKET_VIEW_ASCII: 913 { 914 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure"); 915 } 916 break; 917 918 case DATABUCKET_VIEW_BINARY: 919 { 920 _DataBucketViewBinary(db,filename); 921 } 922 break; 923 924 case DATABUCKET_VIEW_HDF5: 925 { 926 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No HDF5 support"); 927 } 928 break; 929 930 default: 931 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested"); 932 break; 933 } 934 PetscFunctionReturn(0); 935 } 936 937 #undef __FUNCT__ 938 #define __FUNCT__ "DataBucketView_MPI" 939 PetscErrorCode DataBucketView_MPI(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type) 940 { 941 switch (type) { 942 case DATABUCKET_VIEW_STDOUT: 943 { 944 PetscInt f; 945 PetscInt L,buffer,allocated; 946 double memory_usage_total,memory_usage_total_local = 0.0; 947 PetscMPIInt rank; 948 PetscInt ierr; 949 950 ierr = MPI_Comm_rank(comm,&rank); 951 952 DataBucketGetGlobalSizes(comm,db,&L,&buffer,&allocated); 953 954 for( f=0; f<db->nfields; f++ ) { 955 double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6; 956 957 memory_usage_total_local += memory_usage_f; 958 } 959 MPI_Allreduce(&memory_usage_total_local,&memory_usage_total,1,MPI_DOUBLE,MPI_SUM,comm); 960 961 if (rank==0) { 962 PetscPrintf(comm,"DataBucketView(MPI): (\"%s\")\n",filename); 963 PetscPrintf(comm," L = %D \n", L ); 964 PetscPrintf(comm," buffer (max) = %D \n", buffer ); 965 PetscPrintf(comm," allocated = %D \n", allocated ); 966 967 PetscPrintf(comm," nfields registered = %D \n", db->nfields ); 968 for( f=0; f<db->nfields; f++ ) { 969 double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6; 970 971 printf(" [%3d]: field name ==>> %30s : Mem. usage = %1.2e (MB) : rank0\n", f, db->field[f]->name, memory_usage_f ); 972 } 973 974 printf(" Total mem. usage = %1.2e (MB) : collective\n", memory_usage_total ); 975 } 976 977 } 978 break; 979 980 case DATABUCKET_VIEW_ASCII: 981 { 982 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying data structure"); 983 } 984 break; 985 986 case DATABUCKET_VIEW_BINARY: 987 { 988 char *name; 989 PetscMPIInt rank; 990 991 /* create correct extension */ 992 MPI_Comm_rank(comm,&rank); 993 asprintf(&name,"%s_p%1.5d",filename, rank ); 994 995 _DataBucketViewBinary(db,name); 996 997 free(name); 998 } 999 break; 1000 1001 case DATABUCKET_VIEW_HDF5: 1002 { 1003 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for HDF5"); 1004 } 1005 break; 1006 1007 default: 1008 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested"); 1009 break; 1010 } 1011 PetscFunctionReturn(0); 1012 } 1013 1014 #undef __FUNCT__ 1015 #define __FUNCT__ "DataBucketView" 1016 PetscErrorCode DataBucketView(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type) 1017 { 1018 PetscMPIInt nproc; 1019 1020 MPI_Comm_size(comm,&nproc); 1021 if (nproc==1) { 1022 DataBucketView_SEQ(db,filename,type); 1023 } else { 1024 DataBucketView_MPI(comm,db,filename,type); 1025 } 1026 PetscFunctionReturn(0); 1027 } 1028 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "DataBucketDuplicateFields" 1031 PetscErrorCode DataBucketDuplicateFields(DataBucket dbA,DataBucket *dbB) 1032 { 1033 DataBucket db2; 1034 PetscInt f; 1035 1036 DataBucketCreate(&db2); 1037 1038 /* copy contents from dbA into db2 */ 1039 for (f=0; f<dbA->nfields; f++) { 1040 DataField field; 1041 size_t atomic_size; 1042 char *name; 1043 1044 field = dbA->field[f]; 1045 1046 atomic_size = field->atomic_size; 1047 name = field->name; 1048 1049 DataBucketRegisterField(db2,"DataBucketDuplicateFields",name,atomic_size,NULL); 1050 } 1051 DataBucketFinalize(db2); 1052 DataBucketSetInitialSizes(db2,0,1000); 1053 1054 /* set pointer */ 1055 *dbB = db2; 1056 PetscFunctionReturn(0); 1057 } 1058 1059 /* 1060 Insert points from db2 into db1 1061 db1 <<== db2 1062 */ 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "DataBucketInsertValues" 1065 PetscErrorCode DataBucketInsertValues(DataBucket db1,DataBucket db2) 1066 { 1067 PetscInt n_mp_points1,n_mp_points2; 1068 PetscInt n_mp_points1_new,p; 1069 1070 DataBucketGetSizes(db1,&n_mp_points1,0,0); 1071 DataBucketGetSizes(db2,&n_mp_points2,0,0); 1072 1073 n_mp_points1_new = n_mp_points1 + n_mp_points2; 1074 DataBucketSetSizes(db1,n_mp_points1_new,-1); 1075 1076 for (p=0; p<n_mp_points2; p++) { 1077 // db1 <<== db2 // 1078 DataBucketCopyPoint( db2,p, db1,(n_mp_points1 + p) ); 1079 } 1080 PetscFunctionReturn(0); 1081 } 1082 1083 /* helpers for parallel send/recv */ 1084 #undef __FUNCT__ 1085 #define __FUNCT__ "DataBucketCreatePackedArray" 1086 PetscErrorCode DataBucketCreatePackedArray(DataBucket db,size_t *bytes,void **buf) 1087 { 1088 PetscInt f; 1089 size_t sizeof_marker_contents; 1090 void *buffer; 1091 1092 sizeof_marker_contents = 0; 1093 for (f=0; f<db->nfields; f++) { 1094 DataField df = db->field[f]; 1095 1096 sizeof_marker_contents += df->atomic_size; 1097 } 1098 1099 buffer = malloc(sizeof_marker_contents); 1100 memset(buffer,0,sizeof_marker_contents); 1101 1102 if (bytes) { *bytes = sizeof_marker_contents; } 1103 if (buf) { *buf = buffer; } 1104 PetscFunctionReturn(0); 1105 } 1106 1107 #undef __FUNCT__ 1108 #define __FUNCT__ "DataBucketDestroyPackedArray" 1109 PetscErrorCode DataBucketDestroyPackedArray(DataBucket db,void **buf) 1110 { 1111 if (buf) { 1112 free(*buf); 1113 *buf = NULL; 1114 } 1115 PetscFunctionReturn(0); 1116 } 1117 1118 #undef __FUNCT__ 1119 #define __FUNCT__ "DataBucketFillPackedArray" 1120 PetscErrorCode DataBucketFillPackedArray(DataBucket db,const PetscInt index,void *buf) 1121 { 1122 PetscInt f; 1123 void *data,*data_p; 1124 size_t asize,offset; 1125 1126 offset = 0; 1127 for (f=0; f<db->nfields; f++) { 1128 DataField df = db->field[f]; 1129 1130 asize = df->atomic_size; 1131 1132 data = (void*)( df->data ); 1133 data_p = (void*)( (char*)data + index*asize ); 1134 1135 memcpy( (void*)((char*)buf + offset), data_p, asize); 1136 offset = offset + asize; 1137 } 1138 PetscFunctionReturn(0); 1139 } 1140 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "DataBucketInsertPackedArray" 1143 PetscErrorCode DataBucketInsertPackedArray(DataBucket db,const PetscInt idx,void *data) 1144 { 1145 PetscInt f; 1146 void *data_p; 1147 size_t offset; 1148 1149 offset = 0; 1150 for (f=0; f<db->nfields; f++) { 1151 DataField df = db->field[f]; 1152 1153 data_p = (void*)( (char*)data + offset ); 1154 1155 DataFieldInsertPoint(df, idx, (void*)data_p ); 1156 offset = offset + df->atomic_size; 1157 } 1158 PetscFunctionReturn(0); 1159 } 1160