/[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 4705 by jfenwick, Fri Feb 21 02:36:15 2014 UTC revision 4712 by sshaw, Wed Feb 26 04:08:41 2014 UTC
# Line 19  Line 19 
19  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
20  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
21  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
22    #include <ripley/LameAssembler3D.h>
23    #include <ripley/domainhelpers.h>
24  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 43  using esysUtils::FileWriter; Line 45  using esysUtils::FileWriter;
45    
46  namespace ripley {  namespace ripley {
47    
48    int indexOfMax(int a, int b, int c) {
49        if (a > b) {
50            if (c > a) {
51                return 2;
52            }
53            return 0;
54        } else if (b > c) {
55            return 1;
56        }
57        return 2;
58    }
59    
60  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
61               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
62               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
# Line 56  Brick::Brick(int n0, int n1, int n2, dou Line 70  Brick::Brick(int n0, int n1, int n2, dou
70          d2=1;          d2=1;
71      }      }
72      bool warn=false;      bool warn=false;
73      // if number of subdivisions is non-positive, try to subdivide by the same  
74      // ratio as the number of elements      std::vector<int> factors;
75      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
76          warn=true;      int epr[3] = {n0,n1,n2};
77          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
78          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
79          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
80          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
81          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
82              // ratios not the same so leave "smallest" side undivided and try                  continue;
83              // dividing 2 sides only              }
84              if (n0>=n1) {              epr[i] = -1; // can no longer be max
85                  if (n1>=n2) {              if (ranks % d[i] != 0) {
86                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
87                      d2=1;              }
88                  } else {              //remove
89                      d0=d2=0;              ranks /= d[i];
90                      d1=1;          }
91                  }          factorise(factors, ranks);
92              } else {          if (factors.size() != 0) {
93                  if (n0>=n2) {              warn = true;
94                      d0=d1=0;          }
95                      d2=1;      }
96                  } else {      while (factors.size() > 0) {
97                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
98                      d1=d2=0;          int f = factors.back();
99                  }          factors.pop_back();
100              }          d[i] *= f;
101          }          epr[i] /= f;
     }  
     if (d0<=0 && d1<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));  
         d1=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n1) {  
                 d0=0;  
                 d1=1;  
             } else {  
                 d0=1;  
                 d1=0;  
             }  
         }  
     } else if (d0<=0 && d2<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n2) {  
                 d0=0;  
                 d2=1;  
             } else {  
                 d0=1;  
                 d2=0;  
             }  
         }  
     } else if (d1<=0 && d2<=0) {  
         warn=true;  
         d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d1;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n1>n2) {  
                 d1=0;  
                 d2=1;  
             } else {  
                 d1=1;  
                 d2=0;  
             }  
         }  
     }  
     if (d0<=0) {  
         // d1,d2 are preset, determine d0  
         d0=m_mpiInfo->size/(d1*d2);  
     } else if (d1<=0) {  
         // d0,d2 are preset, determine d1  
         d1=m_mpiInfo->size/(d0*d2);  
     } else if (d2<=0) {  
         // d0,d1 are preset, determine d2  
         d2=m_mpiInfo->size/(d0*d1);  
102      }      }
103        d0 = d[0]; d1 = d[1]; d2 = d[2];
104    
105      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
106      // among number of ranks      // among number of ranks
107      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
108          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
109        }
110      if (warn) {      if (warn) {
111          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
112              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
# Line 2867  namespace Line 2829  namespace
2829      {      {
2830          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
2831          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
2832      double total=0;          double total=0;
2833      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
2834      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
2835      {          {
2836          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
2837          {              {
2838          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
2839          {                          {        
2840              arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));                      arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
2841              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];                          total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
2842          }                  }
2843          }              }
2844      }          }
2845      double invtotal=1/total;          double invtotal=1/total;
2846      for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)          for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
2847      {          {
2848          arr[p]*=invtotal;              arr[p]*=invtotal;
2849      }          }
2850      return arr;          return arr;
2851      }      }
2852            
2853      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2893  namespace Line 2855  namespace
2855      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
2856      {      {
2857          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
2858      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
2859      double result=0;          double result=0;
2860      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
2861      {          {
2862          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
2863          {                  {    
2864          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
2865          {                  {
2866              result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];                      result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
2867          }                  }
2868          }              }
2869      }          }
2870      // use this line for "pass-though" (return the centre point value)          // use this line for "pass-though" (return the centre point value)
2871  //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
2872      return result;                return result;      
2873      }      }
2874  }  }
2875    
# Line 2921  escript::Data Brick::randomFill(const es Line 2883  escript::Data Brick::randomFill(const es
2883      int numvals=escript::DataTypes::noValues(shape);      int numvals=escript::DataTypes::noValues(shape);
2884      if (len(filter)>0 && (numvals!=1))      if (len(filter)>0 && (numvals!=1))
2885      {      {
2886      throw RipleyException("Ripley only supports filters for scalar data.");          throw RipleyException("Ripley only supports filters for scalar data.");
2887      }      }
2888      escript::Data res=randomFillWorker(shape, seed, filter);      escript::Data res=randomFillWorker(shape, seed, filter);
2889      if (res.getFunctionSpace()!=what)      if (res.getFunctionSpace()!=what)
2890      {      {
2891      escript::Data r=escript::Data(res, what);          escript::Data r=escript::Data(res, what);
2892      return r;          return r;
2893      }      }
2894      return res;      return res;
2895  }  }
# Line 2976  escript::Data Brick::randomFillWorker(co Line 2938  escript::Data Brick::randomFillWorker(co
2938          throw RipleyException("Brick must be 3D.");          throw RipleyException("Brick must be 3D.");
2939      }      }
2940            
2941      unsigned int radius=0;  // these are only used by gaussian      unsigned int radius=0;  // these are only used by gaussian
2942      double sigma=0.5;      double sigma=0.5;
2943            
2944      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
2945            
2946      if (len(filter)==0)      if (len(filter)==0)
2947      {      {
2948      // nothing special required here yet      // nothing special required here yet
2949      }      }
2950      else if (len(filter)==3)      else if (len(filter)==3)
2951      {      {
2952      boost::python::extract<string> ex(filter[0]);          boost::python::extract<string> ex(filter[0]);
2953      if (!ex.check() || (ex()!="gaussian"))          if (!ex.check() || (ex()!="gaussian"))
2954      {          {
2955          throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
2956      }          }
2957      boost::python::extract<unsigned int> ex1(filter[1]);          boost::python::extract<unsigned int> ex1(filter[1]);
2958      if (!ex1.check())          if (!ex1.check())
2959      {          {
2960          throw RipleyException("Radius of gaussian filter must be a positive integer.");              throw RipleyException("Radius of gaussian filter must be a positive integer.");
2961      }          }
2962      radius=ex1();          radius=ex1();
2963      sigma=0.5;          sigma=0.5;
2964      boost::python::extract<double> ex2(filter[2]);          boost::python::extract<double> ex2(filter[2]);
2965      if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
2966      {          {
2967          throw RipleyException("Sigma must be a postive floating point number.");              throw RipleyException("Sigma must be a postive floating point number.");
2968      }                      }            
2969      }      }
2970      else      else
2971      {      {
# Line 3014  escript::Data Brick::randomFillWorker(co Line 2976  escript::Data Brick::randomFillWorker(co
2976    
2977            
2978      size_t internal[3];      size_t internal[3];
2979      internal[0]=m_NE[0]+1;  // number of points in the internal region      internal[0]=m_NE[0]+1;  // number of points in the internal region
2980      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of
2981      internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of      internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of
2982      size_t ext[3];      size_t ext[3];
2983      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=internal[0]+2*radius;    // includes points we need as input
2984      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=internal[1]+2*radius;    // for smoothing
2985      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=internal[2]+2*radius;    // for smoothing
2986            
2987      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
2988      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
# Line 3045  escript::Data Brick::randomFillWorker(co Line 3007  escript::Data Brick::randomFillWorker(co
3007        
3008      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3009      {      {
3010      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
3011      // will be thrown on all ranks      // will be thrown on all ranks
3012      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
3013    
3014      }      }
3015      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
# Line 3076  cout << "basex=" << basex << " basey=" < Line 3038  cout << "basex=" << basex << " basey=" <
3038    
3039    
3040      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3041      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3042          // there is an overlap of two points both of which need to have "radius" points on either side.          // there is an overlap of two points both of which need to have "radius" points on either side.
3043            
3044      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3045      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3046      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3047            
3048      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3049            
3050      MPI_Request reqs[50];       // a non-tight upper bound on how many we need      MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3051      MPI_Status stats[50];      MPI_Status stats[50];
3052      short rused=0;      short rused=0;
3053            
# Line 3101  cout << "basex=" << basex << " basey=" < Line 3063  cout << "basex=" << basex << " basey=" <
3063      int comserr=0;          int comserr=0;    
3064      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3065      {      {
3066      message& m=incoms[i];          message& m=incoms[i];
3067      comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3068      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3069      }      }
3070    
3071      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3072      {      {
3073      message& m=outcoms[i];          message& m=outcoms[i];
3074      comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3075      }          }    
3076            
3077      if (!comserr)      if (!comserr)
# Line 3119  cout << "basex=" << basex << " basey=" < Line 3081  cout << "basex=" << basex << " basey=" <
3081    
3082      if (comserr)      if (comserr)
3083      {      {
3084      // Yes this is throwing an exception as a result of an MPI error.      // Yes this is throwing an exception as a result of an MPI error.
3085      // and no we don't inform the other ranks that we are doing this.      // and no we don't inform the other ranks that we are doing this.
3086      // however, we have no reason to believe coms work at this point anyway      // however, we have no reason to believe coms work at this point anyway
3087          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3088      }      }
3089            
# Line 3129  cout << "basex=" << basex << " basey=" < Line 3091  cout << "basex=" << basex << " basey=" <
3091            
3092  #endif      #endif    
3093            
3094      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3095      {      {
3096                
3097      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3098      escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3099      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3100      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3101            
3102            
3103      // now we need to copy values over          // now we need to copy values over
3104      for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3105      {          {
3106          for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)    
3107          {              {
3108          for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3109          {                  {
3110              for (unsigned int i=0;i<numvals;++i)                      for (unsigned int i=0;i<numvals;++i)
3111              {                      {
3112                  dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];                          dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3113              }                      }
3114          }                  }
3115          }              }
3116      }            }  
3117      delete[] src;          delete[] src;
3118      return resdat;                return resdat;      
3119      }      }
3120      else        // filter enabled        else        // filter enabled  
3121      {      {
3122            
3123      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3124      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3125      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3126      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3127      double* convolution=get3DGauss(radius, sigma);          double* convolution=get3DGauss(radius, sigma);
3128    
3129      for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3130      {          {
3131          for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)    
3132          {              {
3133          for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3134          {                      {    
3135              dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3136                            
3137          }                  }
3138          }              }
3139      }          }
3140            
3141      delete[] convolution;          delete[] convolution;
3142      delete[] src;          delete[] src;
3143      return resdat;          return resdat;
3144            
3145      }      }
3146  }  }

Legend:
Removed from v.4705  
changed lines
  Added in v.4712

  ViewVC Help
Powered by ViewVC 1.1.26