00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <hdf5.h>
00005 #include "H5Part.h"
00006
00007
00008 herr_t H5PartIOcounter(hid_t group_id,
00009 const char *member_name,
00010 void *operator_data);
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 H5PartFile *H5PartOpenFileParallel(const char *filename,unsigned flags
00032 #ifdef PARALLEL_IO
00033 , MPI_Comm comm
00034 #endif
00035 ){
00036 H5PartFile *f=(H5PartFile*)malloc(sizeof(H5PartFile));
00037 #ifdef PARALLEL_IO
00038 MPI_Info info = MPI_INFO_NULL;
00039 #endif
00040 f->xfer_prop = f->create_prop = f->access_prop=H5P_DEFAULT;
00041 #ifdef PARALLEL_IO
00042 MPI_Comm_size(comm,&f->nprocs);
00043 MPI_Comm_rank(comm,&f->myproc);
00044 f->pnparticles=(long long*)malloc(sizeof(long long)*f->nprocs);
00045
00046 f->access_prop = H5Pcreate(H5P_FILE_ACCESS);
00047 if(H5Pset_fapl_mpio(f->access_prop,comm,info)<0){
00048 fprintf(stderr,"Total failure trying to setup mpi-io for access\n");
00049 exit(0);
00050 }
00051
00052
00053 f->create_prop = H5P_DEFAULT;
00054
00055
00056
00057 f->xfer_prop = H5Pcreate(H5P_DATASET_XFER);
00058 H5Pset_dxpl_mpio(f->xfer_prop,H5FD_MPIO_COLLECTIVE);
00059 f->comm = comm;
00060 #endif
00061 if(flags==H5PART_READ){
00062
00063
00064 f->file = H5Fopen(filename,H5F_ACC_RDONLY,H5P_DEFAULT);
00065 }
00066 else if(flags==H5PART_WRITE){
00067 f->file=H5Fcreate(filename,H5F_ACC_TRUNC,f->create_prop,f->access_prop);
00068 }
00069 else {
00070 f->file=-1;
00071 fprintf(stderr,"H5Open: Invalid file access type %d\n",flags);
00072 }
00073 if(f->file<=0) {
00074 free(f);
00075 fprintf(stderr,"H5Open: File open failed for file [%s]\n",
00076 filename);
00077 exit(0);
00078 return 0;
00079 }
00080 else {
00081 fprintf(stderr,"Proc[%u] Opened file %s val=%d\n",
00082 f->myproc,filename,f->file);
00083 }
00084 f->mode=flags;
00085 f->maxstep=0;
00086 f->timegroup=0;
00087 f->timestep=0;
00088 f->shape=0;
00089 f->diskshape=H5S_ALL;
00090 f->memshape=H5S_ALL;
00091 return f;
00092 }
00093
00094 H5PartFile *H5PartOpenFile(const char *filename,unsigned flags){
00095 H5PartFile *f=(H5PartFile*)malloc(sizeof(H5PartFile));
00096
00097
00098 f->xfer_prop = f->create_prop = f->access_prop=H5P_DEFAULT;
00099 #ifdef PARALLEL_IO
00100 f->pnparticles=0;
00101 f->comm = MPI_COMM_WORLD;
00102 f->nprocs=1;
00103 f->myproc=0;
00104 #endif
00105 if(flags==H5PART_READ){
00106
00107
00108 f->file = H5Fopen(filename,H5F_ACC_RDONLY,H5P_DEFAULT);
00109 if(f->file==-1){
00110 fprintf(stderr,"H5Open: File %s not found\n",filename);
00111 }
00112 }
00113 else if(flags==H5PART_WRITE){
00114 f->file=H5Fcreate(filename,H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
00115 }
00116 else {
00117 f->file=-1;
00118 fprintf(stderr,"H5Open: Invalid file access type %d\n",flags);
00119 }
00120 if(f->file<0) {
00121 free(f);
00122 fprintf(stderr,"H5Open: File open failed for file [%s]\n",
00123 filename);
00124 return 0;
00125 }
00126 f->mode=flags;
00127 f->maxstep=0;
00128 f->timegroup=0;
00129 f->timestep=0;
00130 f->shape=0;
00131 f->diskshape=H5S_ALL;
00132 f->memshape=H5S_ALL;
00133 return f;
00134 }
00135 int H5PartFileIsValid(H5PartFile *f){
00136 if(f->file > 0) return 1;
00137 else return 0;
00138 }
00139
00140
00141
00142 void H5PartCloseFile(H5PartFile *f){
00143 if(f->shape>0) H5Sclose(f->shape); f->shape=0;
00144 if(f->timegroup>0) H5Gclose(f->timegroup); f->timegroup=0;
00145 if(f->diskshape!=H5S_ALL) H5Sclose(f->diskshape);
00146 #ifdef PARALLEL_IO
00147 if(f->pnparticles) free(f->pnparticles);
00148 #endif
00149 if(f->xfer_prop!=H5P_DEFAULT) {
00150 H5Pclose(f->xfer_prop); f->xfer_prop=H5P_DEFAULT;
00151 }
00152 if(f->access_prop!=H5P_DEFAULT) {
00153 H5Pclose(f->access_prop); f->access_prop=H5P_DEFAULT;
00154 }
00155 if(f->create_prop!=H5P_DEFAULT) {
00156 H5Pclose(f->create_prop); f->create_prop=H5P_DEFAULT;
00157 }
00158 if(f->file) H5Fclose(f->file); f->file=0;
00159 free(f);
00160 }
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 void H5PartSetNumParticles(H5PartFile *f,long long nparticles){
00172 hssize_t start[1];
00173 hsize_t stride[1],count[1],total;
00174 register int i,r;
00175
00176 if(f->nparticles==nparticles) return;
00177 if(f->diskshape!=H5S_ALL) H5Sclose(f->diskshape);
00178 if(f->memshape!=H5S_ALL) H5Sclose(f->memshape);
00179 f->memshape=f->diskshape = H5S_ALL;
00180 if(f->shape) H5Sclose(f->shape);
00181 f->nparticles=(hsize_t)nparticles;
00182 #ifndef PARALLEL_IO
00183 f->shape=H5Screate_simple(1,
00184 &(f->nparticles),
00185 NULL);
00186 #else
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 MPI_Allgather(&nparticles,1,MPI_LONG_LONG,f->pnparticles,1,MPI_LONG_LONG,f->comm);
00197
00198 if(f->myproc==0){
00199 puts("AllGather: Particle offsets:\n");
00200 for(i=0;i<f->nprocs;i++)
00201 printf("\tnp=%u\n",f->pnparticles[i]);
00202 }
00203
00204 stride[0]=1;
00205 start[0]=0;
00206 for(i=0;i<f->myproc;i++) start[0]+=f->pnparticles[i];
00207 total=0;
00208 for(i=0;i<f->nprocs;i++) total+=f->pnparticles[i];
00209 f->shape = H5Screate_simple(1,&total,NULL);
00210 f->diskshape = H5Screate_simple(1,&total,NULL);
00211 f->memshape = H5Screate_simple(1,&(f->nparticles),NULL);
00212
00213 if(f->shape<0 || f->memshape<0 || f->diskshape<0){
00214 fprintf(stderr,"Abort: shape construction failed\n");
00215 exit(0);
00216 }
00217 count[0]=nparticles;
00218
00219 r=H5Sselect_hyperslab(f->diskshape,H5S_SELECT_SET,start,stride,count,NULL);
00220 if(r<0){
00221 fprintf(stderr,"Abort: Selection Failed!\n");
00222 exit(0);
00223 }
00224 #endif
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 int H5PartWriteDataFloat64(H5PartFile *f,char *name,double *array){
00237 register int r;
00238 hid_t dataset;
00239
00240
00241 if(f->timegroup<=0) {fprintf(stderr,"Eeeerrrrroooorrrr!!!!\n");}
00242 dataset = H5Dcreate(f->timegroup,name,H5T_NATIVE_DOUBLE,f->shape,H5P_DEFAULT);
00243 if(dataset<0)
00244 fprintf(stderr,"Dataset open failed for name=[%s] step=%u!\n",
00245 name,f->timestep);
00246 r=H5Dwrite(dataset,H5T_NATIVE_DOUBLE,f->memshape,f->diskshape,H5P_DEFAULT,array);
00247 H5Dclose(dataset);
00248 return r;
00249 }
00250 int H5PartWriteDataInt64(H5PartFile *f,char *name,long long *array){
00251 register int r;
00252 hid_t dataset;
00253
00254
00255 if(f->timegroup<=0) {fprintf(stderr,"Eeeerrrrroooorrrr!!!!\n");}
00256 dataset = H5Dcreate(f->timegroup,name,H5T_NATIVE_INT64,f->shape,H5P_DEFAULT);
00257 if(dataset<0){
00258 fprintf(stderr,"Dataset open failed for name=[%s] step=%u!\n",
00259 name,f->timestep);
00260 exit(0);
00261 }
00262 r=H5Dwrite(dataset,H5T_NATIVE_INT64,f->memshape,f->diskshape,H5P_DEFAULT,array);
00263 if(r<0) {
00264 fprintf(stderr,"Attempt to write dataset failed!\n");
00265 exit(0);
00266 }
00267 H5Dclose(dataset);
00268 return r;
00269 }
00270 int H5PartWriteFileAttribString(H5PartFile *f,char *name,
00271 char *attrib){
00272 return H5PartWriteFileAttrib(f,name,H5T_NATIVE_CHAR,attrib,strlen(attrib)+1);
00273 }
00274 int H5PartWriteStepAttribString(H5PartFile *f,char *name,
00275 char *attrib){
00276 return H5PartWriteStepAttrib(f,name,H5T_NATIVE_CHAR,attrib,strlen(attrib)+1);
00277 }
00278 int H5PartWriteStepAttrib(H5PartFile *f,char *name,
00279 hid_t type,void *value,int nelem){
00280 register int r;
00281 hid_t attrib;
00282 hid_t space;
00283 hsize_t len;
00284
00285
00286 if(f->timegroup<=0) {fprintf(stderr,"Eeeerrrrroooorrrr!!!!\n");}
00287 len = nelem;
00288 space = H5Screate_simple(1,&len,NULL);
00289 attrib = H5Acreate(f->timegroup,name,type,space,H5P_DEFAULT);
00290 r=H5Awrite(attrib,type,value);
00291 H5Aclose(attrib);
00292 H5Sclose(space);
00293 return r;
00294 }
00295
00296 int H5PartWriteAttrib(H5PartFile *f,char *name,
00297 hid_t type,void *value,int nelem){
00298 return H5PartWriteStepAttrib(f,name,type,value,nelem);
00299 }
00300
00301
00302 int H5PartWriteFileAttrib(H5PartFile *f,char *name,
00303 hid_t type,void *value,int nelem){
00304 register int r;
00305 hid_t attrib;
00306 hid_t space;
00307 hsize_t len;
00308
00309
00310 if(f->file<=0) {fprintf(stderr,"Eeeerrrrroooorrrr!!!!\n");}
00311 len = nelem;
00312 space = H5Screate_simple(1,&len,NULL);
00313 attrib = H5Acreate(f->file,name,type,space,H5P_DEFAULT);
00314 r=H5Awrite(attrib,type,value);
00315 H5Aclose(attrib);
00316 H5Sclose(space);
00317 return r;
00318 }
00319 herr_t H5PartAttribcounter(hid_t group_id,
00320 const char *member_name,
00321 void *operator_data){
00322 int *count = (int*)operator_data;
00323 (*count)++;
00324 return 0;
00325 }
00326 int H5PartGetNumStepAttribs(H5PartFile *f){
00327 return H5Aget_num_attrs(f->timegroup);
00328 }
00329 int H5PartGetNumFileAttribs(H5PartFile *f){
00330
00331
00332 int nattribs;
00333 hid_t rootgroup=H5Gopen(f->file,"/");;
00334 nattribs=H5Aget_num_attrs(rootgroup);
00335
00336
00337 H5Gclose(rootgroup);
00338 return nattribs;
00339 }
00340
00341 hid_t H5PartNormType(hid_t type){
00342 H5T_class_t tclass = H5Tget_class(type);
00343 int size = H5Tget_size(type);
00344 switch(tclass){
00345 case H5T_INTEGER:
00346 if(size==8) return H5T_NATIVE_INT64;
00347 else if(size==1) return H5T_NATIVE_CHAR;
00348 else fprintf(stderr,"Unknown type %u. Cannot infer normalized type\n",type);
00349 break;
00350 case H5T_FLOAT:
00351 return H5T_NATIVE_DOUBLE;
00352 break;
00353 }
00354 return -1;
00355 }
00356
00357 void H5PartGetStepAttribInfo(H5PartFile *f,int idx,
00358 char *name,size_t maxname,hid_t *type,int *nelem){
00359 hid_t attrib=H5Aopen_idx(f->timegroup,idx);
00360 hid_t mytype=H5Aget_type(attrib);
00361 hid_t space = H5Aget_space(attrib);
00362 if(nelem)
00363 *nelem=H5Sget_simple_extent_npoints(space);
00364 if(name)
00365 H5Aget_name(attrib,maxname,name);
00366 if(type)
00367 *type=H5PartNormType(mytype);
00368
00369
00370
00371
00372 H5Sclose(space);
00373 H5Tclose(mytype);
00374 H5Aclose(attrib);
00375 }
00376 void H5PartGetFileAttribInfo(H5PartFile *f,int idx,
00377 char *name,size_t maxname,hid_t *type,int *nelem){
00378 hid_t attrib=H5Aopen_idx(f->file,idx);
00379 hid_t mytype=H5Aget_type(attrib);
00380 hid_t space = H5Aget_space(attrib);
00381 if(nelem)
00382 *nelem=H5Sget_simple_extent_npoints(space);
00383 if(name)
00384 H5Aget_name(attrib,maxname,name);
00385 if(type)
00386 *type=H5PartNormType(mytype);
00387
00388
00389
00390
00391 H5Sclose(space);
00392 H5Tclose(mytype);
00393 H5Aclose(attrib);
00394 }
00395
00396 void H5PartReadStepAttrib(H5PartFile *f,char *name,void *data){
00397
00398 hid_t attrib=H5Aopen_name(f->timegroup,name);
00399 hid_t mytype=H5Aget_type(attrib);
00400 hid_t space = H5Aget_space(attrib);
00401 hsize_t nelem=H5Sget_simple_extent_npoints(space);
00402 hid_t type=H5PartNormType(mytype);
00403
00404
00405
00406
00407 H5Aread(attrib,type,data);
00408 H5Sclose(space);
00409 H5Tclose(mytype);
00410 H5Aclose(attrib);
00411 }
00412
00413 void H5PartReadAttrib(H5PartFile *f,char *name,void *data){
00414
00415 H5PartReadStepAttrib(f,name,data);
00416 }
00417
00418
00419
00420
00421 void H5PartReadFileAttrib(H5PartFile *f,char *name,void *data){
00422 hid_t attrib=H5Aopen_name(f->file,name);
00423 hid_t mytype=H5Aget_type(attrib);
00424 hid_t space = H5Aget_space(attrib);
00425 hsize_t nelem=H5Sget_simple_extent_npoints(space);
00426 hid_t type=H5PartNormType(mytype);
00427
00428
00429
00430
00431 H5Aread(attrib,type,data);
00432 H5Sclose(space);
00433 H5Tclose(mytype);
00434 H5Aclose(attrib);
00435 }
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 void H5PartSetStep(H5PartFile *f,int step){
00450 char name[128];
00451
00452 f->timestep=step;
00453 if(f->timegroup>0) {
00454 H5Gclose(f->timegroup);
00455 }
00456 f->timegroup=-1;
00457 sprintf(name,"Particles#%u",f->timestep);
00458
00459
00460 if(f->mode==H5PART_READ) {
00461
00462 f->timegroup=H5Gopen(f->file,name);
00463 }
00464 if(f->timegroup<=0){
00465
00466
00467 f->timegroup = H5Gcreate(f->file,name,0);
00468 if(f->timegroup<0)
00469 fprintf(stderr,"p[%u] timegroup creation failed!\n",
00470 f->myproc);
00471 }
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481 herr_t H5PartIOcounter(hid_t group_id,
00482 const char *member_name,
00483 void *operator_data){
00484 int *count = (int*)operator_data;
00485 H5G_stat_t objinfo;
00486
00487 if(!strncmp(member_name,"Particles",9)) (*count)++;
00488 return 0;
00489 }
00490 herr_t H5PartDScounter(hid_t group_id,
00491 const char *member_name,
00492 void *operator_data){
00493 int *count = (int*)operator_data;
00494 H5G_stat_t objinfo;
00495
00496 if(H5Gget_objinfo(group_id,
00497 member_name,
00498 1 ,
00499 &objinfo)<0) {
00500 return 0;
00501 }
00502 if(objinfo.type==H5G_DATASET) (*count)++;
00503 return 0;
00504 }
00505
00506 typedef struct H5IO_getname_t {
00507 int index,count;
00508 char *name;
00509 }H5IO_getname_t;
00510
00511 herr_t H5IOgetname(hid_t group_id,
00512 const char *member_name,
00513 void *operator_data){
00514 H5IO_getname_t *getn = (H5IO_getname_t*)operator_data;
00515
00516 H5G_stat_t objinfo;
00517
00518 if(H5Gget_objinfo(group_id,
00519 member_name,
00520 1 ,
00521 &objinfo)<0) return 0;
00522
00523 if(objinfo.type!=H5G_DATASET)
00524 return 0;
00525 if(getn->index==getn->count){
00526 strcpy(getn->name,member_name);
00527 return 1;
00528 }
00529 getn->count++;
00530 return 0;
00531 }
00532
00533
00534
00535
00536
00537
00538
00539 int H5PartGetNumSteps(H5PartFile *f){
00540
00541 int count=0;
00542 int idx=0;
00543 while(H5Giterate(f->file,
00544 "/",
00545 &idx,
00546 H5PartIOcounter,
00547 &count)<0){}
00548 return count;
00549 }
00550 int H5PartGetNumDatasets(H5PartFile *f){
00551 int count=0;
00552 int idx=0;
00553 char name[128];
00554 sprintf(name,"Particles#%u",f->timestep);
00555 while(H5Giterate(f->file,
00556 name,
00557 &idx,
00558 H5PartDScounter,
00559 &count)<0){}
00560 return count;
00561 }
00562 int H5PartGetDatasetName(H5PartFile *f,
00563 int _index,
00564 char *name,
00565 int maxlen){
00566 H5IO_getname_t getn;
00567 int idx=_index;
00568 char gname[128];
00569 sprintf(gname,"Particles#%u",f->timestep);
00570 getn.index=_index; getn.name=name; getn.count=_index;
00571 while(H5Giterate(f->file,
00572 gname,
00573 &idx,
00574 H5IOgetname,
00575 &getn)<0){}
00576 return 1;
00577 }
00578
00579
00580
00581
00582
00583
00584
00585 long long H5PartGetNumParticles(H5PartFile *f){
00586 hid_t space,dataset;
00587 hsize_t nparticles;
00588 if(f->timegroup<=0) H5PartSetStep(f,f->timestep);
00589
00590 dataset=H5Dopen(f->timegroup,"x");
00591 space = H5Dget_space(dataset);
00592 nparticles= H5Sget_simple_extent_npoints(space);
00593 H5Sclose(space);
00594 H5Dclose(dataset);
00595 return (long long)nparticles;
00596 }
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 int H5PartReadDataFloat64(H5PartFile *f,char *name,double *array){
00607 hid_t space,dataset,datatype;
00608 if(!f->timegroup) H5PartSetStep(f,f->timestep);
00609 dataset=H5Dopen(f->timegroup,name);
00610 space = H5Dget_space(dataset);
00611
00612 H5Dread(dataset,
00613 H5T_NATIVE_DOUBLE,
00614
00615 H5S_ALL,
00616 H5S_ALL,
00617 H5P_DEFAULT,
00618 array);
00619 H5Sclose(space);
00620 H5Dclose(dataset);
00621 return 1;
00622 }
00623
00624 int H5PartReadDataInt64(H5PartFile *f,char *name,long long *array){
00625 hid_t space,dataset,datatype;
00626 if(!f->timegroup) H5PartSetStep(f,f->timestep);
00627 dataset=H5Dopen(f->timegroup,name);
00628 space = H5Dget_space(dataset);
00629
00630 H5Dread(dataset,
00631 H5T_NATIVE_INT64,
00632
00633 H5S_ALL,
00634 H5S_ALL,
00635 H5P_DEFAULT,
00636 array);
00637 H5Sclose(space);
00638 H5Dclose(dataset);
00639 return 1;
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 int H5PartReadParticleStep(H5PartFile *f,
00652 int step,
00653 double *x,double *y,double *z,
00654 double *px,double *py,double *pz,
00655 long long *id){
00656 H5PartSetStep(f,step);
00657
00658 H5PartReadDataFloat64(f,"x",x);
00659 H5PartReadDataFloat64(f,"y",y);
00660 H5PartReadDataFloat64(f,"z",z);
00661 H5PartReadDataFloat64(f,"px",px);
00662 H5PartReadDataFloat64(f,"py",py);
00663 H5PartReadDataFloat64(f,"pz",pz);
00664 H5PartReadDataInt64(f,"id",id);
00665 return 1;
00666 }