# Diff of /branches/diaplayground/ripley/src/Brick.cpp

revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC revision 4277 by caltinay, Wed Mar 6 01:30:41 2013 UTC
# Line 1  Line 1
1
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
# Line 18  extern "C" { Line 18  extern "C" {
18  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
19  }  }
20
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24
25  #if USE_SILO  #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 54  Brick::Brick(int n0, int n1, int n2, dou Line 58  Brick::Brick(int n0, int n1, int n2, dou
58      if (d0<=0 && d1<=0 && d2<=0) {      if (d0<=0 && d1<=0 && d2<=0) {
59          warn=true;          warn=true;
60          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));
61          d1=(int)(d0*n1/(float)n0);          d0=max(1, d0);
62            d1=max(1, (int)(d0*n1/(float)n0));
63          d2=m_mpiInfo->size/(d0*d1);          d2=m_mpiInfo->size/(d0*d1);
64          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
65              // ratios not the same so leave "smallest" side undivided and try              // ratios not the same so leave "smallest" side undivided and try
# Line 80  Brick::Brick(int n0, int n1, int n2, dou Line 85  Brick::Brick(int n0, int n1, int n2, dou
85      }      }
86      if (d0<=0 && d1<=0) {      if (d0<=0 && d1<=0) {
87          warn=true;          warn=true;
88          d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));          d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
89          d1=m_mpiInfo->size/d0;          d1=m_mpiInfo->size/d0;
90          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
91              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 94  Brick::Brick(int n0, int n1, int n2, dou Line 99  Brick::Brick(int n0, int n1, int n2, dou
99          }          }
100      } else if (d0<=0 && d2<=0) {      } else if (d0<=0 && d2<=0) {
101          warn=true;          warn=true;
102          d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1)));          d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));
103          d2=m_mpiInfo->size/d0;          d2=m_mpiInfo->size/d0;
104          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
105              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 108  Brick::Brick(int n0, int n1, int n2, dou Line 113  Brick::Brick(int n0, int n1, int n2, dou
113          }          }
114      } else if (d1<=0 && d2<=0) {      } else if (d1<=0 && d2<=0) {
115          warn=true;          warn=true;
116          d1=int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1)));          d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));
117          d2=m_mpiInfo->size/d1;          d2=m_mpiInfo->size/d1;
118          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
119              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 242  bool Brick::operator==(const AbstractDom Line 247  bool Brick::operator==(const AbstractDom
247
249                             const vector<int>& first,                             const vector<int>& first,
250                             const vector<int>& numValues) const                             const vector<int>& numValues,
251                               const vector<int>& multiplier) const
252  {  {
253      // check destination function space      // check destination function space
254      int myN0, myN1, myN2;      int myN0, myN1, myN2;
264      } else      } else
265          throw RipleyException("readBinaryGrid(): invalid function space for output data object");          throw RipleyException("readBinaryGrid(): invalid function space for output data object");
266
267        if (first.size() != 3)
268            throw RipleyException("readBinaryGrid(): argument 'first' must have 3 entries");
269
270        if (numValues.size() != 3)
271            throw RipleyException("readBinaryGrid(): argument 'numValues' must have 3 entries");
272
273        if (multiplier.size() != 3)
274            throw RipleyException("readBinaryGrid(): argument 'multiplier' must have 3 entries");
275        for (size_t i=0; i<multiplier.size(); i++)
276            if (multiplier[i]<1)
277                throw RipleyException("readBinaryGrid(): all multipliers must be positive");
278
279      // check file existence and size      // check file existence and size
280      ifstream f(filename.c_str(), ifstream::binary);      ifstream f(filename.c_str(), ifstream::binary);
281      if (f.fail()) {      if (f.fail()) {
291      }      }
292
293      // check if this rank contributes anything      // check if this rank contributes anything
294      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset0 ||
295              first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||              first[1] >= m_offset1+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset1 ||
296              first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {              first[2] >= m_offset2+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset2) {
297          f.close();          f.close();
298          return;          return;
299      }      }
308      const int idx0 = max(0, m_offset0-first[0]);      const int idx0 = max(0, m_offset0-first[0]);
309      const int idx1 = max(0, m_offset1-first[1]);      const int idx1 = max(0, m_offset1-first[1]);
310      const int idx2 = max(0, m_offset2-first[2]);      const int idx2 = max(0, m_offset2-first[2]);
311      // number of values to write      // number of values to read
312      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(numValues[0]-idx0, myN0-first0);
313      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(numValues[1]-idx1, myN1-first1);
314      const int num2 = min(numValues[2]-idx2, myN2-first2);      const int num2 = min(numValues[2]-idx2, myN2-first2);
322              const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);              const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
323              f.seekg(fileofs*sizeof(float));              f.seekg(fileofs*sizeof(float));
325
326              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
327                  double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0+(first2+z)*myN0*myN1);                  const int baseIndex = first0+x*multiplier[0]
328                  for (index_t c=0; c<numComp; c++) {                                          +(first1+y*multiplier[1])*myN0
329                      for (index_t q=0; q<dpp; q++) {                                          +(first2+z*multiplier[2])*myN0*myN1;
330                          *dest++ = static_cast<double>(values[x*numComp+c]);                  for (index_t m2=0; m2<multiplier[2]; m2++) {
331                        for (index_t m1=0; m1<multiplier[1]; m1++) {
332                            for (index_t m0=0; m0<multiplier[0]; m0++) {
333                                const int dataIndex = baseIndex+m0
334                                               +m1*myN0
335                                               +m2*myN0*myN1;
336                                double* dest = out.getSampleDataRW(dataIndex);
337                                for (index_t c=0; c<numComp; c++) {
338                                    if (!isnan(values[x*numComp+c])) {
339                                        for (index_t q=0; q<dpp; q++) {
340                                            *dest++ = static_cast<double>(values[x*numComp+c]);
341                                        }
342                                    }
343                                }
344                            }
345                      }                      }
346                  }                  }
347              }              }
351      f.close();      f.close();
352  }  }
353
354    void Brick::readNcGrid(escript::Data& out, string filename, string varname,
355                const vector<int>& first, const vector<int>& numValues,
356                const vector<int>& multiplier) const
357    {
358    #ifdef USE_NETCDF
359        // check destination function space
360        int myN0, myN1, myN2;
361        if (out.getFunctionSpace().getTypeCode() == Nodes) {
362            myN0 = m_N0;
363            myN1 = m_N1;
364            myN2 = m_N2;
365        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
366                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
367            myN0 = m_NE0;
368            myN1 = m_NE1;
369            myN2 = m_NE2;
370        } else
371            throw RipleyException("readNcGrid(): invalid function space for output data object");
372
373        if (first.size() != 3)
374            throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
375
376        if (numValues.size() != 3)
377            throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
378
379        if (multiplier.size() != 3)
380            throw RipleyException("readNcGrid(): argument 'multiplier' must have 3 entries");
381        for (size_t i=0; i<multiplier.size(); i++)
382            if (multiplier[i]<1)
383                throw RipleyException("readNcGrid(): all multipliers must be positive");
384
385        // check file existence and size
387        if (!f.is_valid())
388            throw RipleyException("readNcGrid(): cannot open file");
389
390        NcVar* var = f.get_var(varname.c_str());
391        if (!var)
392            throw RipleyException("readNcGrid(): invalid variable name");
393
394        // TODO: rank>0 data support
395        const int numComp = out.getDataPointSize();
396        if (numComp > 1)
397            throw RipleyException("readNcGrid(): only scalar data supported");
398
399        const int dims = var->num_dims();
400        const long *edges = var->edges();
401
402        // is this a slice of the data object (dims!=3)?
403        // note the expected ordering of edges (as in numpy: z,y,x)
404        if ( (dims==3 && (numValues[2] > edges[0] || numValues[1] > edges[1]
405                          || numValues[0] > edges[2]))
406                || (dims==2 && numValues[2]>1)
407                || (dims==1 && (numValues[2]>1 || numValues[1]>1)) ) {
408            throw RipleyException("readNcGrid(): not enough data in file");
409        }
410
411        // check if this rank contributes anything
412        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset0 ||
413                first[1] >= m_offset1+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset1 ||
414                first[2] >= m_offset2+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset2) {
415            return;
416        }
417
418        // now determine how much this rank has to write
419
420        // first coordinates in data object to write to
421        const int first0 = max(0, first[0]-m_offset0);
422        const int first1 = max(0, first[1]-m_offset1);
423        const int first2 = max(0, first[2]-m_offset2);
424        // indices to first value in file
425        const int idx0 = max(0, m_offset0-first[0]);
426        const int idx1 = max(0, m_offset1-first[1]);
427        const int idx2 = max(0, m_offset2-first[2]);
428        // number of values to read
429        const int num0 = min(numValues[0]-idx0, myN0-first0);
430        const int num1 = min(numValues[1]-idx1, myN1-first1);
431        const int num2 = min(numValues[2]-idx2, myN2-first2);
432
433        vector<double> values(num0*num1*num2);
434        if (dims==3) {
435            var->set_cur(idx2, idx1, idx0);
436            var->get(&values[0], num2, num1, num0);
437        } else if (dims==2) {
438            var->set_cur(idx1, idx0);
439            var->get(&values[0], num1, num0);
440        } else {
441            var->set_cur(idx0);
442            var->get(&values[0], num0);
443        }
444
445        const int dpp = out.getNumDataPointsPerSample();
446        out.requireWrite();
447
448        for (index_t z=0; z<num2; z++) {
449            for (index_t y=0; y<num1; y++) {
450    #pragma omp parallel for
451                for (index_t x=0; x<num0; x++) {
452                    const int baseIndex = first0+x*multiplier[0]
453                                            +(first1+y*multiplier[1])*myN0
454                                            +(first2+z*multiplier[2])*myN0*myN1;
455                    const int srcIndex=z*num1*num0+y*num0+x;
456                    if (!isnan(values[srcIndex])) {
457                        for (index_t m2=0; m2<multiplier[2]; m2++) {
458                            for (index_t m1=0; m1<multiplier[1]; m1++) {
459                                for (index_t m0=0; m0<multiplier[0]; m0++) {
460                                    const int dataIndex = baseIndex+m0
461                                                   +m1*myN0
462                                                   +m2*myN0*myN1;
463                                    double* dest = out.getSampleDataRW(dataIndex);
464                                    for (index_t q=0; q<dpp; q++) {
465                                        *dest++ = values[srcIndex];
466                                    }
467                                }
468                            }
469                        }
470                    }
471                }
472            }
473        }
474    #else
475        throw RipleyException("readNcGrid(): not compiled with netCDF support");
476    #endif
477    }
478
479  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
480  {  {
481  #if USE_SILO  #if USE_SILO
# Line 1758  void Brick::dofToNodes(escript::Data& ou Line 1916  void Brick::dofToNodes(escript::Data& ou
1916                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1917          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1918      }      }
1919        Paso_Coupler_free(coupler);
1920  }  }
1921
1922  //private  //private

Legend:
 Removed from v.3981 changed lines Added in v.4277