fermiqcd: /Users/mdipierro/fermiqcd/development/Libraries/mdp_field_save

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00014 template<class T>
00015 bool mdp_field<T>::save_vtk(string filename, 
00016                             int t,
00017                             int component,
00018                             int processIO,
00019                             bool ASCII) {
00020 
00021   filename=next_to_latest_file(filename);
00022   string filename_tmp=filename+".tmp";
00023 
00024   int max_buffer_size=1024;
00025   int timeslice;
00026   mdp_int header_size=0;
00027   mdp_int psize=field_components*Tsize;
00028   mdp_int idx_gl, nvol_gl=lattice().nvol_gl, k;
00029   double mytime=mpi.time();
00030   header.reset();
00031   if(lattice().ndim<3 && lattice().ndim>4) error("mdp_field::save_vtk only works for ndim=3 and 4");
00032   if(ME==processIO) {
00033     mdp_int *buffer_size=new mdp_int[Nproc];
00034     mdp_int *buffer_ptr =new mdp_int[Nproc];
00035     mdp_array<T,3> large_buffer(Nproc,max_buffer_size,field_components);
00036     T *short_buffer=new T[field_components];
00037     int process;
00038     for(process=0; process<Nproc; process++) buffer_ptr[process]=0;
00039     cout << "Saving file " << filename 
00040          << " from process " << processIO 
00041          << " (buffer = " << max_buffer_size << " sites)" << '\n'; 
00042     fflush(stdout);
00043     FILE *fp=0;    
00044 
00045     float fval;
00046     char tmp[1024];
00047     char header[1024];
00048     int skip_bytes=0;
00049         
00050     int space_volume=lattice().size()/lattice().size(0);
00051     int LZ=lattice().size(1);
00052     int LY=lattice().size(2);
00053     int LX=lattice().size(3);    
00054     if(lattice().ndim==3) {
00055       space_volume=lattice().size();
00056       LZ=lattice().size(0);
00057       LY=lattice().size(1);
00058       LX=lattice().size(2);
00059       t=-1;
00060     }
00061     int fc,fc2;
00062     if(component==-1) fc2=fc=field_components; else fc2=fc=component;
00063 
00064     sprintf(header,
00065             "# vtk DataFile Version 2.0\n"
00066             "%s\n"
00067             "%s\n"
00068             "DATASET STRUCTURED_POINTS\n"
00069             "DIMENSIONS %i %i %i\n"
00070             "ORIGIN     0   0   0\n"
00071             "SPACING    1   1   1\n"
00072             "POINT_DATA %i",
00073             filename.c_str(),
00074             ((ASCII)?"ASCII":"BINARY"),
00075             LX,LY,LZ,LX*LY*LZ);
00076 
00077     fp=fopen(filename_tmp.c_str(),"wb");
00078     fwrite(header,sizeof(char),strlen(header),fp);
00079    
00080     for(idx_gl=0; idx_gl<nvol_gl; idx_gl++) {
00081       process=where_global(idx_gl);
00082       if((process!=NOWHERE) && (process!=processIO)) {
00083         if(buffer_ptr[process]==0) {
00084           mpi.get(buffer_size[process], process);
00085           mpi.get(&(large_buffer(process,0,0)), 
00086                   buffer_size[process]*field_components, process);
00087         }
00088         for(k=0; k<field_components; k++)
00089           short_buffer[k]=large_buffer(process,buffer_ptr[process],k);
00090         buffer_ptr[process]++;
00091         if(buffer_ptr[process]==buffer_size[process]) buffer_ptr[process]=0; 
00092       }
00093       if(process==processIO) {
00094         for(k=0; k<field_components; k++)
00095           short_buffer[k]=*(m+lattice().local(idx_gl)*field_components+k);
00096       }
00097       if(process!=NOWHERE) {
00098         timeslice=idx_gl/space_volume;
00099         if(idx_gl % space_volume==0 && (t<0 || timeslice==t)) {   
00100           sprintf(header,"\nSCALARS scalars_t%i float\nLOOKUP_TABLE default\n",timeslice);
00101           fwrite(header,sizeof(char),strlen(header),fp);
00102         }
00103         if(t<0 || timeslice==t || lattice().ndim==3) {    
00104           for(fc=0; fc<field_components; fc++) if(component==-1 || fc==component) {
00105             fval=(float) short_buffer[fc];
00106             if(!ASCII) {
00107               switch_endianess_byte4(fval);
00108               if(fwrite(&fval,sizeof(fval),1,fp)!=1) 
00109                 error("probably out of disk space");
00110             } else {
00111               sprintf(tmp,"%e ",fval);
00112               if(fwrite(tmp,strlen(tmp), 1,fp)!=1)
00113                 error("probably out of disk space");
00114             }
00115           }
00116         }
00117       }
00118     }
00119     if(fp) {
00120       fclose(fp);
00121       remove(filename.c_str());
00122       rename(filename_tmp.c_str(),filename.c_str());
00123     }
00124     delete[] buffer_size;
00125     delete[] buffer_ptr;
00126     delete[] short_buffer;
00127   } else {
00128     int  process;
00129     mdp_int buffer_size=0, idx, idx_gl;
00130     mdp_int *local_index=new mdp_int[max_buffer_size];
00131     mdp_array<T,2> local_buffer(max_buffer_size,field_components);
00132     mdp_request request;
00133     for(idx_gl=0; idx_gl<nvol_gl; idx_gl++) {
00134       process=where_global(idx_gl);
00135       if(process==ME) {
00136         local_index[buffer_size]=lattice().local(idx_gl);
00137         buffer_size++;
00138       }
00139       if((buffer_size==max_buffer_size) || 
00140          ((idx_gl==nvol_gl-1) && (buffer_size>0))) {
00141         for(idx=0; idx<buffer_size; idx++)
00142           for(k=0; k<field_components; k++)
00143             local_buffer(idx,k)=*(m+local_index[idx]*field_components+k);
00144         mpi.put(buffer_size, processIO, request);
00145         mpi.wait(request);
00146         mpi.put(&(local_buffer(0,0)), buffer_size*field_components, 
00147                 processIO, request);
00148         mpi.wait(request);
00149         buffer_size=0;
00150       }
00151     }
00152     delete[] local_index;
00153   }
00154   if(ME==0 && mdp_shutup==false) {
00155     printf("... Saving time: %f (sec)\n", mpi.time()-mytime);
00156     fflush(stdout);
00157   }
00158   return true;
00159 }
00160 
00161 mdp_field<float>& cumulate_field(mdp_field<float>& field, string filename) {
00162   static map<string,int> counter;
00163   static map<string,mdp_field<float>*> fields;
00164   mdp_site p(field.lattice());
00165   int k=0;
00166   if(counter.find(filename)==counter.end()) {
00167     fields[filename]=new mdp_field<float>(field.lattice(),field.size_per_site());
00168     forallsites(p) for(int i=0; i<field.size_per_site(); i++) (*fields[filename])(p,i)=field(p,i);
00169     counter[filename]=1;
00170   } else {
00171     k=counter[filename];
00172     forallsites(p) for(int i=0; i<field.size_per_site(); i++)
00173       (*fields[filename])(p,i)=(k*(*fields[filename])(p,i)+field(p,i))/(k+1);
00174     counter[filename]++;
00175   }
00176   return (*fields[filename]);
00177 }
00178 
00179 /*
00180 void save_vtk(mdp_field<float> &field,
00181               string filename, 
00182               int t=-1,
00183               int component=-1,
00184               int processIO=0,
00185               bool ASCII=false) {
00186   static map<string,int> counter;
00187   static map<string,mdp_field<float>*> fields;
00188   mdp_site p(field.lattice());
00189   int k=0;
00190   if(counter.find(filename)==counter.end()) {
00191     fields[filename]=new mdp_field<float>(field.lattice(),field.size_per_site());
00192     forallsites(p) for(int i=0; i<field.size_per_site(); i++) (*fields[filename])(p,i)=field(p,i);
00193     counter[filename]=1;    
00194   } else {
00195     k=counter[filename];
00196     forallsites(p) for(int i=0; i<field.size_per_site(); i++) 
00197       (*fields[filename])(p,i)=(k*(*fields[filename])(p,i)+field(p,i))/(k+1);
00198     counter[filename]++;
00199   }
00200   if(filename[filename.size()-1]!='*') {
00201     filename=filename.substr(0,filename.size()-4);
00202     field.save_vtk(filename+".vtk",t,component,processIO,ASCII);
00203     fields[filename]->save_vtk(filename+"_average.vtk",t,component,processIO,ASCII);
00204   } else {
00205     string filename2=filename.substr(0,filename.size()-1);
00206     field.save_vtk(filename2+"_"+tostring(k)+".vtk",t,component,processIO,ASCII);
00207     fields[filename]->save_vtk(filename2+"_average_"+tostring(k)+".vtk",t,component,processIO,ASCII);
00208   }
00209 }
00210 */
/Users/mdipierro/fermiqcd/development/Libraries/mdp_field_save_vtk.h
