/[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 4650 by jfenwick, Wed Feb 5 04:16:01 2014 UTC revision 4705 by jfenwick, Fri Feb 21 02:36:15 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
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)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
# Line 786  const int* Brick::borrowSampleReferenceI Line 787  const int* Brick::borrowSampleReferenceI
787          case FaceElements:          case FaceElements:
788          case ReducedFaceElements:          case ReducedFaceElements:
789              return &m_faceId[0];              return &m_faceId[0];
790            case Points:
791                return &m_diracPointNodeIDs[0];
792          default:          default:
793              break;              break;
794      }      }
# Line 2872  namespace Line 2875  namespace
2875          {          {
2876          for (int x=-r;x<=r;++x)          for (int x=-r;x<=r;++x)
2877          {                  {        
2878              arr[(x+r)+(y+r)*(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));
2879              total+=arr[(x+r)+(y+r)*(r*2+1)];              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
2880          }          }
2881          }          }
2882      }      }
2883      double invtotal=1/total;      double invtotal=1/total;
2884      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)      for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
2885      {      {
2886          arr[p]*=invtotal;          arr[p]*=invtotal;
2887      }      }
# Line 2902  namespace Line 2905  namespace
2905          }          }
2906          }          }
2907      }      }
2908          return result;            // use this line for "pass-though" (return the centre point value)
2909    //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];
2910        return result;      
2911      }      }
2912  }  }
2913    
2914    /* This is a wrapper for filtered (and non-filtered) randoms
2915     * For detailed doco see randomFillWorker
2916    */
2917    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
2918           const escript::FunctionSpace& what,
2919           long seed, const boost::python::tuple& filter) const
2920    {
2921        int numvals=escript::DataTypes::noValues(shape);
2922        if (len(filter)>0 && (numvals!=1))
2923        {
2924        throw RipleyException("Ripley only supports filters for scalar data.");
2925        }
2926        escript::Data res=randomFillWorker(shape, seed, filter);
2927        if (res.getFunctionSpace()!=what)
2928        {
2929        escript::Data r=escript::Data(res, what);
2930        return r;
2931        }
2932        return res;
2933    }
2934    
2935  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
2936  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
# Line 2944  inset=2*radius+1 Line 2969  inset=2*radius+1
2969  This is to ensure that values at distance `radius` from the shared/overlapped element  This is to ensure that values at distance `radius` from the shared/overlapped element
2970  that ripley has.  that ripley has.
2971  */  */
2972  escript::Data Brick::randomFill(long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
2973  {  {
2974      if (m_numDim!=3)      if (m_numDim!=3)
2975      {      {
2976          throw RipleyException("Brick must be 3D.");          throw RipleyException("Brick must be 3D.");
2977      }      }
2978      if (len(filter)!=3) {      
2979          throw RipleyException("Unsupported random filter");      unsigned int radius=0;  // these are only used by gaussian
2980      }      double sigma=0.5;
2981      boost::python::extract<string> ex(filter[0]);      
2982      if (!ex.check() || (ex()!="gaussian"))      unsigned int numvals=escript::DataTypes::noValues(shape);
2983        
2984        if (len(filter)==0)
2985      {      {
2986          throw RipleyException("Unsupported random filter");      // nothing special required here yet
2987      }      }
2988      boost::python::extract<unsigned int> ex1(filter[1]);      else if (len(filter)==3)
     if (!ex1.check())  
2989      {      {
2990          throw RipleyException("Radius of gaussian filter must be a positive integer.");      boost::python::extract<string> ex(filter[0]);
2991        if (!ex.check() || (ex()!="gaussian"))
2992        {
2993            throw RipleyException("Unsupported random filter for Brick.");
2994        }
2995        boost::python::extract<unsigned int> ex1(filter[1]);
2996        if (!ex1.check())
2997        {
2998            throw RipleyException("Radius of gaussian filter must be a positive integer.");
2999        }
3000        radius=ex1();
3001        sigma=0.5;
3002        boost::python::extract<double> ex2(filter[2]);
3003        if (!ex2.check() || (sigma=ex2())<=0)
3004        {
3005            throw RipleyException("Sigma must be a postive floating point number.");
3006        }            
3007      }      }
3008      unsigned int radius=ex1();      else
     double sigma=0.5;  
     boost::python::extract<double> ex2(filter[2]);  
     if (!ex2.check() || (sigma=ex2())<=0)  
3009      {      {
3010          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Unsupported random filter");
3011      }          }
3012        
3013        
3014    
3015            
3016      size_t internal[3];      size_t internal[3];
3017      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
# Line 2983  escript::Data Brick::randomFill(long see Line 3025  escript::Data Brick::randomFill(long see
3025      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3026      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3027    
3028      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3029      {      {
3030          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3031      }      }
3032            if (2*radius>=internal[1]-4)
3033        {
3034            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3035        }
3036        if (2*radius>=internal[2]-4)
3037        {
3038            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3039        }    
3040        
3041      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3042      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3043            
     
3044  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3045          
3046        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3047        {
3048        // since the dimensions are equal for all ranks, this exception
3049        // will be thrown on all ranks
3050        throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
3051    
3052        }
3053      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3054      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3055      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3056    #endif    
3057    
3058    /*    
3059        // if we wanted to test a repeating pattern
3060        size_t basex=0;
3061        size_t basey=0;
3062        size_t basez=0;
3063    #ifdef ESYS_MPI    
3064        basex=X*m_gNE[0]/m_NX[0];
3065        basey=Y*m_gNE[1]/m_NX[1];
3066        basez=Z*m_gNE[2]/m_NX[2];
3067            
3068    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3069        
3070    #endif    
3071        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3072    */
3073    
3074    #ifdef ESYS_MPI
3075    
3076    
3077    
3078      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);
3079      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3080            // there is an overlap of two points both of which need to have "radius" points on either side.
3081            
3082      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
3083      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3084      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3085            
3086      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3087            
3088      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
3089      MPI_Status stats[50];      MPI_Status stats[50];
# Line 3019  escript::Data Brick::randomFill(long see Line 3096  escript::Data Brick::randomFill(long see
3096      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3097            
3098            
3099      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
       
3100            
3101      int comserr=0;          int comserr=0;    
3102      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
# Line 3030  escript::Data Brick::randomFill(long see Line 3106  escript::Data Brick::randomFill(long see
3106      block.setUsed(m.destbuffid);      block.setUsed(m.destbuffid);
3107      }      }
3108    
3109      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3110      {      {
3111      message& m=incoms[i];      message& m=outcoms[i];
3112      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++));
3113      }          }    
3114            
# Line 3050  escript::Data Brick::randomFill(long see Line 3126  escript::Data Brick::randomFill(long see
3126      }      }
3127            
3128      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3129        
3130  #endif      #endif    
3131      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3132      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
     // don't need to check for exwrite because we just made it  
     escript::DataVector& dv=resdat.getExpandedVectorReference();  
     double* convolution=get3DGauss(radius, sigma);  
     for (size_t z=0;z<(internal[2]);++z)  
3133      {      {
3134      for (size_t y=0;y<(internal[1]);++y)            
3135        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3136        escript::Data resdat(0, shape, fs , true);
3137        // don't need to check for exwrite because we just made it
3138        escript::DataVector& dv=resdat.getExpandedVectorReference();
3139        
3140        
3141        // now we need to copy values over
3142        for (size_t z=0;z<(internal[2]);++z)
3143      {      {
3144          for (size_t x=0;x<(internal[0]);++x)          for (size_t y=0;y<(internal[1]);++y)    
3145          {              {
3146          dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);          for (size_t x=0;x<(internal[0]);++x)
3147                    {
3148                for (unsigned int i=0;i<numvals;++i)
3149                {
3150                    dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3151                }
3152            }
3153            }
3154        }  
3155        delete[] src;
3156        return resdat;      
3157        }
3158        else        // filter enabled  
3159        {
3160        
3161        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3162        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3163        // don't need to check for exwrite because we just made it
3164        escript::DataVector& dv=resdat.getExpandedVectorReference();
3165        double* convolution=get3DGauss(radius, sigma);
3166    
3167        for (size_t z=0;z<(internal[2]);++z)
3168        {
3169            for (size_t y=0;y<(internal[1]);++y)    
3170            {
3171            for (size_t x=0;x<(internal[0]);++x)
3172            {    
3173                dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3174                
3175            }
3176          }          }
3177      }      }
3178        
3179        delete[] convolution;
3180        delete[] src;
3181        return resdat;
3182        
3183      }      }
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3184  }  }
3185    
3186    
# Line 3081  escript::Data Brick::randomFill(long see Line 3191  escript::Data Brick::randomFill(long see
3191  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3192      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3193      //is the found element even owned by this rank      //is the found element even owned by this rank
3194        // (inside owned or shared elements but will map to an owned element)
3195      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3196          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3197                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3198            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3199                    + m_dx[dim]/2.;
3200            if (min > coords[dim] || max < coords[dim]) {
3201              return NOT_MINE;              return NOT_MINE;
3202          }          }
3203      }      }
# Line 3102  int Brick::findNode(const double *coords Line 3216  int Brick::findNode(const double *coords
3216          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3217      }      }
3218      //find the closest node      //find the closest node
3219      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3220          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3221          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3222              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3223              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3224                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3225                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3226                  if (total < minDist) {                  if (total < minDist) {
3227                      closest = INDEX3(ex+dy, ey+dy, ez+dz, m_NE[0]+1, m_NE[1]+1);                      closest = INDEX3(ex+dy-m_offset[0], ey+dy-m_offset[1],
3228                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3229                      minDist = total;                      minDist = total;
3230                  }                  }
3231              }              }
3232          }          }
3233      }      }
3234        if (closest == NOT_MINE) {
3235            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3236        }
3237      return closest;      return closest;
3238  }  }
3239    

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

  ViewVC Help
Powered by ViewVC 1.1.26