/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3971 by caltinay, Wed Sep 19 02:55:35 2012 UTC revision 4174 by caltinay, Wed Jan 30 03:21:27 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  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <ripley/Brick.h>  #include <ripley/Brick.h>
17  extern "C" {  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 52  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 78  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 92  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 106  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 302  void Brick::readBinaryGrid(escript::Data Line 309  void Brick::readBinaryGrid(escript::Data
309              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]);
310              f.seekg(fileofs*sizeof(float));              f.seekg(fileofs*sizeof(float));
311              f.read((char*)&values[0], num0*numComp*sizeof(float));              f.read((char*)&values[0], num0*numComp*sizeof(float));
312                const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1;
313    
314              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
315                  double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0+(first2+z)*myN0*myN1);                  double* dest = out.getSampleDataRW(dataIndex+x);
316                  for (index_t c=0; c<numComp; c++) {                  for (index_t c=0; c<numComp; c++) {
317                      for (index_t q=0; q<dpp; q++) {                      if (!isnan(values[x*numComp+c])) {
318                          *dest++ = static_cast<double>(values[x*numComp+c]);                          for (index_t q=0; q<dpp; q++) {
319                                *dest++ = static_cast<double>(values[x*numComp+c]);
320                            }
321                      }                      }
322                  }                  }
323              }              }
# Line 316  void Brick::readBinaryGrid(escript::Data Line 327  void Brick::readBinaryGrid(escript::Data
327      f.close();      f.close();
328  }  }
329    
330    void Brick::readNcGrid(escript::Data& out, string filename, string varname,
331                const vector<int>& first, const vector<int>& numValues) const
332    {
333    #ifdef USE_NETCDF
334        // check destination function space
335        int myN0, myN1, myN2;
336        if (out.getFunctionSpace().getTypeCode() == Nodes) {
337            myN0 = m_N0;
338            myN1 = m_N1;
339            myN2 = m_N2;
340        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
341                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
342            myN0 = m_NE0;
343            myN1 = m_NE1;
344            myN2 = m_NE2;
345        } else
346            throw RipleyException("readNcGrid(): invalid function space for output data object");
347    
348        if (first.size() != 3)
349            throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
350    
351        if (numValues.size() != 3)
352            throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
353    
354        // check file existence and size
355        NcFile f(filename.c_str(), NcFile::ReadOnly);
356        if (!f.is_valid())
357            throw RipleyException("readNcGrid(): cannot open file");
358    
359        NcVar* var = f.get_var(varname.c_str());
360        if (!var)
361            throw RipleyException("readNcGrid(): invalid variable name");
362    
363        // TODO: rank>0 data support
364        const int numComp = out.getDataPointSize();
365        if (numComp > 1)
366            throw RipleyException("readNcGrid(): only scalar data supported");
367    
368        const int dims = var->num_dims();
369        const long *edges = var->edges();
370    
371        // is this a slice of the data object (dims!=3)?
372        // note the expected ordering of edges (as in numpy: z,y,x)
373        if ( (dims==3 && (numValues[2] > edges[0] || numValues[1] > edges[1]
374                          || numValues[0] > edges[2]))
375                || (dims==2 && numValues[2]>1)
376                || (dims==1 && (numValues[2]>1 || numValues[1]>1)) ) {
377            throw RipleyException("readNcGrid(): not enough data in file");
378        }
379    
380        // check if this rank contributes anything
381        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
382                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
383                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
384            return;
385        }
386    
387        // now determine how much this rank has to write
388    
389        // first coordinates in data object to write to
390        const int first0 = max(0, first[0]-m_offset0);
391        const int first1 = max(0, first[1]-m_offset1);
392        const int first2 = max(0, first[2]-m_offset2);
393        // indices to first value in file
394        const int idx0 = max(0, m_offset0-first[0]);
395        const int idx1 = max(0, m_offset1-first[1]);
396        const int idx2 = max(0, m_offset2-first[2]);
397        // number of values to write
398        const int num0 = min(numValues[0]-idx0, myN0-first0);
399        const int num1 = min(numValues[1]-idx1, myN1-first1);
400        const int num2 = min(numValues[2]-idx2, myN2-first2);
401    
402        vector<double> values(num0*num1*num2);
403        if (dims==3) {
404            var->set_cur(idx2, idx1, idx0);
405            var->get(&values[0], num2, num1, num0);
406        } else if (dims==2) {
407            var->set_cur(idx1, idx0);
408            var->get(&values[0], num1, num0);
409        } else {
410            var->set_cur(idx0);
411            var->get(&values[0], num0);
412        }
413    
414        const int dpp = out.getNumDataPointsPerSample();
415        out.requireWrite();
416    
417        for (index_t z=0; z<num2; z++) {
418            for (index_t y=0; y<num1; y++) {
419    #pragma omp parallel for
420                for (index_t x=0; x<num0; x++) {
421                    const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1+x;
422                    const int srcIndex=z*num1*num0+y*num0+x;
423                    if (!isnan(values[srcIndex])) {
424                        double* dest = out.getSampleDataRW(dataIndex);
425                        for (index_t q=0; q<dpp; q++) {
426                            *dest++ = values[srcIndex];
427                        }
428                    }
429                }
430            }
431        }
432    #else
433        throw RipleyException("readNcGrid(): not compiled with netCDF support");
434    #endif
435    }
436    
437  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
438  {  {
439  #if USE_SILO  #if USE_SILO
# Line 1756  void Brick::dofToNodes(escript::Data& ou Line 1874  void Brick::dofToNodes(escript::Data& ou
1874                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1875          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1876      }      }
1877        Paso_Coupler_free(coupler);
1878  }  }
1879    
1880  //private  //private

Legend:
Removed from v.3971  
changed lines
  Added in v.4174

  ViewVC Help
Powered by ViewVC 1.1.26