00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "Utility/FieldBlock.h"
00028 #include "Utility/IpplInfo.h"
00029 #include "Utility/PAssert.h"
00030 #include "Field/BrickExpression.h"
00031 #include "Field/Field.h"
00032 #include "FieldLayout/FieldLayout.h"
00033 #include "Field/LField.h"
00034 #include "Message/Message.h"
00035 #include "Message/Communicate.h"
00036 #include "Profile/Profiler.h"
00037
00038 #include <string.h>
00039 #include <stdio.h>
00040
00041 #ifdef IPPL_NETCDF
00042 #include <netcdf.h>
00043 #endif
00044
00045
00046
00047
00048
00049
00050
00051 template<class T, unsigned Dim, class Mesh, class Centering >
00052 FieldBlock<T,Dim,Mesh,Centering>::FieldBlock(char* fname,
00053 FieldLayout<Dim>& layout,
00054 unsigned numFields) :
00055 Layout(layout), NumFields(numFields)
00056 {
00057 TAU_TYPE_STRING(taustr, CT(*this) + " void (char *, " + CT(layout)
00058 + ", unsigned )" );
00059 TAU_PROFILE("FieldBlock::FieldBlock()", taustr,
00060 TAU_UTILITY | TAU_FIELD | TAU_IO);
00061
00062 strcpy(FName, fname);
00063 NumRecords = 0;
00064
00065
00066 #ifdef IPPL_NETCDF
00067 int Parent = 0;
00068 if( Ippl::Comm->myNode() == Parent ) {
00069 int ncid;
00070 int rcode;
00071
00072
00073 ncid = nccreate(FName, NC_CLOBBER);
00074
00075
00076
00077 int fillmode = ncsetfill(ncid,NC_NOFILL);
00078
00079
00080
00081 const NDIndex<Dim> &domain = Layout.getDomain();
00082
00083
00084 int dimids[4];
00085 if ( Dim==1 ) {
00086 dimids[0] = ncdimdef(ncid, "record", NC_UNLIMITED);
00087 dimids[1] = ncdimdef(ncid, "nx", domain[0].length());
00088 }
00089 if ( Dim==2 ) {
00090 dimids[0] = ncdimdef(ncid, "record", NC_UNLIMITED);
00091 dimids[1] = ncdimdef(ncid, "nx", domain[0].length());
00092 dimids[2] = ncdimdef(ncid, "ny", domain[1].length());
00093 }
00094 if ( Dim==3 ) {
00095 dimids[0] = ncdimdef(ncid, "record", NC_UNLIMITED);
00096 dimids[1] = ncdimdef(ncid, "nx", domain[0].length());
00097 dimids[2] = ncdimdef(ncid, "ny", domain[1].length());
00098 dimids[3] = ncdimdef(ncid, "nz", domain[2].length());
00099 }
00100 if ( Dim>3 ) {
00101 ERRORMSG("FieldBlock: can't write more than 3 dimensions" << endl);
00102 }
00103 for( int i = 0 ; i < NumFields ; i++ ) {
00104 char varName[256];
00105 sprintf(varName, "var%d", i);
00106 ncvardef(ncid, varName, NC_DOUBLE, Dim+1, dimids);
00107 }
00108
00109 int errRet = ncendef(ncid);
00110 rcode = ncclose(ncid);
00111 if ( rcode != 0) {
00112 ERRORMSG("FieldBlock: ncclose() error, rcode=" << rcode << endl);
00113 }
00114 }
00115 #endif // IPPL_NETCDF
00116 }
00117
00118
00119
00120 template<class T, unsigned Dim, class Mesh, class Centering >
00121 FieldBlock<T,Dim,Mesh,Centering>::
00122 FieldBlock(char* fname, FieldLayout<Dim>& layout)
00123 : Layout(layout)
00124 {
00125 TAU_TYPE_STRING(taustr, CT(*this) + " void (char *, " + CT(layout) + " )" );
00126 TAU_PROFILE("FieldBlock::FieldBlock()", taustr,
00127 TAU_UTILITY | TAU_FIELD | TAU_IO);
00128
00129 strcpy(FName, fname);
00130
00131 #ifdef IPPL_NETCDF
00132 int Parent = 0;
00133 if( Ippl::Comm->myNode() == Parent ) {
00134 int i;
00135 int ncid, ndims, nvars;
00136 nc_type datatype;
00137 int rcode;
00138 long netSize[Dim+1];
00139
00140
00141 ncid = ncopen(FName, NC_NOWRITE);
00142
00143
00144 ncinquire(ncid, &ndims, &nvars, 0, 0);
00145 if( ndims != Dim+1)
00146 ERRORMSG("FieldBlock: bad number of dims on read of " << FName << endl);
00147 NumFields = nvars;
00148 for( i = 0; i < Dim+1 ; i++ ) {
00149 ncdiminq(ncid, i, 0, &netSize[i]);
00150 }
00151 for( i = 0; i < Dim ; i++ ) {
00152 if( netSize[i+1] != Layout.getDomain()[i].length() ) {
00153 ERRORMSG("FieldBlock: encountered non-conforming size on axis ");
00154 ERRORMSG(i << endl);
00155 }
00156 }
00157
00158 for( i = 0; i < NumFields ; i++ ) {
00159 ncvarinq(ncid, i, 0, &datatype, 0, 0, 0);
00160 if( datatype != NC_DOUBLE)
00161 ERRORMSG("FieldBlock: file must contain double precion data" << endl);
00162 }
00163 long numRecords;
00164 ncdiminq(ncid, 0, 0, &numRecords);
00165 NumRecords = numRecords;
00166 rcode = ncclose(ncid);
00167 if ( rcode != 0) {
00168 ERRORMSG("FieldBlock: ncclose() error, rcode=" << rcode << endl);
00169 }
00170 }
00171 #endif // IPPL_NETCDF
00172 }
00173
00174 template< class T, unsigned Dim, class Mesh, class Centering >
00175 void FieldBlock<T,Dim,Mesh,Centering>::write(Field<T,Dim,Mesh,Centering>& f,
00176 unsigned varID,
00177 unsigned record)
00178 {
00179 TAU_TYPE_STRING(taustr, "void (" + CT(f) + ", unsigned, unsigned )" );
00180 TAU_PROFILE("FieldBlock::write()", taustr,
00181 TAU_UTILITY | TAU_FIELD | TAU_IO);
00182
00183 #ifdef IPPL_NETCDF
00184 Inform msg("FieldBlock::write");
00185 msg.setPrintNode();
00186
00187 if( varID >=NumFields ) {
00188 ERRORMSG(varID << " is a bad variable ID in FieldBlock::write " << endl);
00189 return;
00190 }
00191 int ncid;
00192 int Parent = 0;
00193
00194
00195 if( Ippl::Comm->myNode() == Parent ) {
00196 ncid = ncopen(FName, NC_WRITE);
00197
00198
00199
00200 int fillmode = ncsetfill(ncid,NC_NOFILL);
00201 }
00202 Ippl::Comm->barrier();
00203
00204 long startIndex[Dim+1];
00205 long countIndex[Dim+1];
00206 int rcode;
00207
00208 int tag = Ippl::Comm->next_tag( FB_WRITE_TAG, FB_TAG_CYCLE );
00209 typedef LField<T,Dim>::iterator LFI;
00210
00211
00212 Field<T,Dim,Mesh,Centering>::iterator_if local;
00213
00214 for (local = f.begin_if(); local != f.end_if(); ++local) {
00215
00216 LField<T,Dim>& l = *(*local).second.get();
00217 NDIndex<Dim>& lo = (NDIndex<Dim>&) l.getOwned();
00218 NDIndex<Dim>& la = (NDIndex<Dim>&) l.getAllocated();
00219
00220
00221 Message *mess = new Message();
00222 ::putMessage(*mess, lo);
00223 LFI msgval = l.begin();
00224 msgval.TryCompress();
00225 ::putMessage(*mess, msgval);
00226
00227
00228
00229
00230 Ippl::Comm->send(mess, Parent, tag);
00231 }
00232
00233
00234
00235
00236 if( Ippl::Comm->myNode() == Parent ) {
00237
00238
00239 int numVnodes = Layout.size_iv() + Layout.size_rdv();
00240
00241
00242 for (int remaining = numVnodes; remaining>0; --remaining) {
00243
00244 int any_node = COMM_ANY_NODE;
00245 Message *mess = Ippl::Comm->receive_block(any_node, tag);
00246 PAssert(mess != 0);
00247
00248
00249
00250
00251
00252 NDIndex<Dim> localBlock;
00253 LFI rhs;
00254 ::getMessage(*mess, localBlock);
00255 ::getMessage(*mess, rhs);
00256
00257
00258 int size = 1;
00259 for( int i = 0 ; i < Dim ; i++ ) {
00260 startIndex[Dim - i] = localBlock[Dim - 1 - i].first();
00261 countIndex[Dim - i] = localBlock[Dim - 1 - i].length();
00262 size *= countIndex[Dim - i];
00263 }
00264
00265
00266 startIndex[0] = record;
00267 countIndex[0] = 1;
00268 double* buffer = new double[size];
00269
00270 int icount = 0;
00271 int n0,n1,n2,i0,i1,i2;
00272 if (rhs.IsCompressed()) {
00273 for (i0=0; i0<size; ++i0) buffer[icount++] = *rhs;
00274
00275 } else {
00276 switch(Dim) {
00277 case 1:
00278 n0 = rhs.size(0);
00279 for (i0=0; i0<n0; ++i0)
00280 buffer[icount++] = rhs.offset(i0);
00281 break;
00282 case 2:
00283 n0 = rhs.size(0);
00284 n1 = rhs.size(1);
00285 for (i0=0; i0<n0; ++i0)
00286 for (i1=0; i1<n1; ++i1)
00287 buffer[icount++] = rhs.offset(i0,i1);
00288 break;
00289 case 3:
00290 n0 = rhs.size(0);
00291 n1 = rhs.size(1);
00292 n2 = rhs.size(2);
00293 for (i0=0; i0<n0; ++i0)
00294 for (i1=0; i1<n1; ++i1)
00295 for (i2=0; i2<n2; ++i2)
00296 buffer[icount++] = rhs.offset(i0,i1,i2);
00297 break;
00298 default:
00299 ERRORMSG("FieldBlock: bad Dimension in Field::write()" << endl);
00300 break;
00301 }
00302 }
00303
00304 rcode = ncvarput(ncid, varID, startIndex,countIndex, buffer);
00305
00306 if ( rcode != 0) {
00307 ERRORMSG("FieldBlock: ncvarput() error, rcode=" << rcode << endl);
00308 }
00309 delete [] buffer;
00310 delete mess;
00311 }
00312 rcode = ncclose(ncid);
00313 if ( rcode != 0) {
00314 ERRORMSG("FieldBlock: ncclose() error, rcode=" << rcode << endl);
00315 }
00316 }
00317 NumRecords = NumRecords > record ? NumRecords : record;
00318 #endif // IPPL_NETCDF
00319 }
00320
00321 template< class T, unsigned Dim, class Mesh, class Centering >
00322 void FieldBlock<T,Dim,Mesh,Centering>::read(Field<T,Dim,Mesh,Centering>& f,
00323 unsigned varID,
00324 unsigned record)
00325 {
00326 TAU_TYPE_STRING(taustr, "void (" + CT(f) + ", unsigned, unsigned )" );
00327 TAU_PROFILE("FieldBlock::read()", taustr,
00328 TAU_UTILITY | TAU_FIELD | TAU_IO);
00329 #ifdef IPPL_NETCDF
00330 Inform msg("FieldBlock::read");
00331 int i;
00332 int ncid;
00333 int Parent = 0;
00334
00335
00336 if( Ippl::Comm->myNode() == Parent ) {
00337 ncid = ncopen(FName, NC_NOWRITE);
00338 }
00339
00340 if( record >= NumRecords )
00341 ERRORMSG("invalid record on FieldBlock::read() " << endl);
00342
00343
00344 long startIndex[Dim+1];
00345 long countIndex[Dim+1];
00346 int rcode;
00347
00348 int tag = Ippl::Comm->next_tag( FB_READ_TAG, FB_TAG_CYCLE );
00349 typedef LField<T,Dim>::iterator LFI;
00350
00351
00352
00353
00354 if( Ippl::Comm->myNode() == Parent ) {
00355
00356
00357
00358 FieldLayout<Dim>::iterator_iv localVnode;
00359 for (localVnode = Layout.begin_iv() ;
00360 localVnode != Layout.end_iv(); ++localVnode) {
00361
00362 Vnode<Dim>& vn = *(*localVnode).second;
00363 NDIndex<Dim>& lo = (NDIndex<Dim>&) vn.getDomain();
00364
00365 int size = 1;
00366 for( i = 0 ; i < Dim ; i++ ) {
00367 startIndex[Dim - i] = lo[Dim - 1 - i].first();
00368 countIndex[Dim - i] = lo[Dim - 1 - i].length();
00369 size *= countIndex[Dim - i];
00370 }
00371
00372 startIndex[0] = record;
00373 countIndex[0] = 1;
00374 double* buffer = new double[size];
00375
00376 rcode = ncvarget(ncid, varID, startIndex, countIndex, buffer);
00377 if ( rcode != 0) {
00378 ERRORMSG("FieldBlock: ncvarget() error, rcode=" << rcode << endl);
00379 }
00380
00381 Message *mess = new Message();
00382 ::putMessage(*mess, lo);
00383 LFI msgval(buffer, lo, lo);
00384 ::putMessage(*mess, msgval);
00385
00386
00387
00388 Ippl::Comm->send(mess, vn.getNode(), tag);
00389 delete [] buffer;
00390 }
00391
00392
00393 FieldLayout<Dim>::iterator_dv remoteVnode;
00394 for (remoteVnode = Layout.begin_rdv() ;
00395 remoteVnode != Layout.end_rdv(); ++remoteVnode) {
00396
00397 NDIndex<Dim>& lo = (NDIndex<Dim>&) (*remoteVnode).first;
00398 Vnode<Dim>& vn = *(*remoteVnode).second;
00399
00400 int size = 1;
00401 for( i = 0 ; i < Dim ; i++ ) {
00402 startIndex[Dim - i] = lo[Dim - 1 - i].first();
00403 countIndex[Dim - i] = lo[Dim - 1 - i].length();
00404 size *= countIndex[Dim - i];
00405 }
00406
00407 startIndex[0] = record;
00408 countIndex[0] = 1;
00409 double* buffer = new double[size];
00410
00411 rcode = ncvarget(ncid, varID, startIndex,countIndex, buffer);
00412 if ( rcode != 0) {
00413 ERRORMSG("FieldBlock: ncvarget() error, rcode=" << rcode << endl);
00414 }
00415
00416
00417 Message *mess = new Message();
00418 ::putMessage(*mess, lo);
00419 LFI msgval(buffer, lo, lo);
00420 ::putMessage(*mess, msgval);
00421
00422
00423 Ippl::Comm->send(mess, vn.getNode(), tag);
00424 delete [] buffer;
00425 }
00426 rcode = ncclose(ncid);
00427 if ( rcode != 0) {
00428 ERRORMSG("FieldBlock: ncclose() error, rcode=" << rcode << endl);
00429 }
00430 }
00431
00432
00433 int numLocalVnodes = Layout.size_iv();
00434 for (int remaining = numLocalVnodes; remaining>0; --remaining) {
00435
00436 int any_node = COMM_ANY_NODE;
00437 Message *mess = Ippl::Comm->receive_block(any_node, tag);
00438 PAssert(mess != 0);
00439
00440
00441 NDIndex<Dim> localBlock;
00442 LFI rhs;
00443 ::getMessage(*mess, localBlock);
00444 ::getMessage(*mess, rhs);
00445
00446
00447 LField<T,Dim>* myLField;
00448 bool flag = 1;
00449
00450 Field<T,Dim,Mesh,Centering>::iterator_if local;
00451 for (local = f.begin_if(); local != f.end_if(); ++local) {
00452 myLField = (*local).second.get();
00453 const NDIndex<Dim>& lo = myLField->getOwned();
00454 if( lo == localBlock ) {
00455 flag = 0;
00456 break;
00457 }
00458 }
00459 if(flag)
00460 ERRORMSG("FieldBlock::read: did not match the local NDIndex" << endl);
00461
00462
00463 LFI lhs = myLField->begin();
00464 int n0,n1,n2,i0,i1,i2;
00465 if( lhs.IsCompressed() && rhs.IsCompressed() ) {
00466 *lhs = *rhs;
00467 }
00468 if( !lhs.IsCompressed() && rhs.IsCompressed() ) {
00469 for (i0=0; i0 < localBlock.size(); ++i0) {
00470 *lhs = *rhs;
00471 ++lhs;
00472 }
00473 }
00474 if( !rhs.IsCompressed() ) {
00475 if( lhs.IsCompressed() ) {
00476 myLField->Uncompress();
00477 lhs = myLField->begin();
00478 }
00479 switch(Dim) {
00480 case 1:
00481 n0 = rhs.size(0);
00482 for (i0=0; i0<n0; ++i0) {
00483 lhs.offset(i0) = *rhs;
00484 ++rhs;
00485 }
00486 break;
00487 case 2:
00488 n0 = rhs.size(0);
00489 n1 = rhs.size(1);
00490 for (i0=0; i0<n0; ++i0)
00491 for (i1=0; i1<n1; ++i1) {
00492 lhs.offset(i0,i1) = *rhs;
00493 ++rhs;
00494 }
00495 break;
00496 case 3:
00497 n0 = rhs.size(0);
00498 n1 = rhs.size(1);
00499 n2 = rhs.size(2);
00500 for (i0=0; i0<n0; ++i0)
00501 for (i1=0; i1<n1; ++i1)
00502 for (i2=0; i2<n2; ++i2) {
00503 lhs.offset(i0,i1,i2) = *rhs;
00504 ++rhs;
00505 }
00506 break;
00507 default:
00508 ERRORMSG("FieldBlock: bad Dimension in write()" << endl);
00509 break;
00510 }
00511 }
00512
00513
00514
00515 delete mess;
00516 }
00517
00518 #endif // IPPL_NETCDF
00519 }
00520
00521
00522
00523
00524
00525