/[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 4703 by jfenwick, Fri Feb 21 00:55:31 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    cerr <<     arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)] << " ";      
2881          }          }
2882    cerr << endl;      
2883          }          }
2884      }      }
2885      double invtotal=1/total;      double invtotal=1/total;
2886      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)
2887      {      {
2888          arr[p]*=invtotal;          arr[p]*=invtotal;
2889      }      }
# Line 2902  namespace Line 2907  namespace
2907          }          }
2908          }          }
2909      }      }
2910          return result;            // use this line for "pass-though" (return the centre point value)
2911    //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];
2912        return result;      
2913      }      }
2914  }  }
2915    
2916    /* This is a wrapper for filtered (and non-filtered) randoms
2917     * For detailed doco see randomFillWorker
2918    */
2919    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
2920           const escript::FunctionSpace& what,
2921           long seed, const boost::python::tuple& filter) const
2922    {
2923        int numvals=escript::DataTypes::noValues(shape);
2924        if (len(filter)>0 && (numvals!=1))
2925        {
2926        throw RipleyException("Ripley only supports filters for scalar data.");
2927        }
2928        escript::Data res=randomFillWorker(shape, seed, filter);
2929        if (res.getFunctionSpace()!=what)
2930        {
2931        escript::Data r=escript::Data(res, what);
2932        return r;
2933        }
2934        return res;
2935    }
2936    
2937  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
2938  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 2971  inset=2*radius+1
2971  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
2972  that ripley has.  that ripley has.
2973  */  */
2974  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
2975  {  {
2976      if (m_numDim!=3)      if (m_numDim!=3)
2977      {      {
2978          throw RipleyException("Brick must be 3D.");          throw RipleyException("Brick must be 3D.");
2979      }      }
2980      if (len(filter)!=3) {      
2981          throw RipleyException("Unsupported random filter");      unsigned int radius=0;  // these are only used by gaussian
2982      }      double sigma=0.5;
2983      boost::python::extract<string> ex(filter[0]);      
2984      if (!ex.check() || (ex()!="gaussian"))      unsigned int numvals=escript::DataTypes::noValues(shape);
2985        
2986        if (len(filter)==0)
2987      {      {
2988          throw RipleyException("Unsupported random filter");      // nothing special required here yet
2989      }      }
2990      boost::python::extract<unsigned int> ex1(filter[1]);      else if (len(filter)==3)
     if (!ex1.check())  
2991      {      {
2992          throw RipleyException("Radius of gaussian filter must be a positive integer.");      boost::python::extract<string> ex(filter[0]);
2993        if (!ex.check() || (ex()!="gaussian"))
2994        {
2995            throw RipleyException("Unsupported random filter for Brick.");
2996        }
2997        boost::python::extract<unsigned int> ex1(filter[1]);
2998        if (!ex1.check())
2999        {
3000            throw RipleyException("Radius of gaussian filter must be a positive integer.");
3001        }
3002        radius=ex1();
3003        sigma=0.5;
3004        boost::python::extract<double> ex2(filter[2]);
3005        if (!ex2.check() || (sigma=ex2())<=0)
3006        {
3007            throw RipleyException("Sigma must be a postive floating point number.");
3008        }            
3009      }      }
3010      unsigned int radius=ex1();      else
     double sigma=0.5;  
     boost::python::extract<double> ex2(filter[2]);  
     if (!ex2.check() || (sigma=ex2())<=0)  
3011      {      {
3012          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Unsupported random filter");
3013      }          }
3014        
3015        
3016    
3017            
3018      size_t internal[3];      size_t internal[3];
3019      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 3027  escript::Data Brick::randomFill(long see
3027      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3028      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3029    
3030      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3031      {      {
3032          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");
3033      }      }
3034            if (2*radius>=internal[1]-4)
3035        {
3036            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3037        }
3038        if (2*radius>=internal[2]-4)
3039        {
3040            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3041        }    
3042        
3043      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3044      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3045            
     
3046  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3047          
3048        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3049        {
3050        // since the dimensions are equal for all ranks, this exception
3051        // will be thrown on all ranks
3052        throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
3053    
3054        }
3055      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3056      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];
3057      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3058    #endif    
3059    
3060        
3061        // if we wanted to test a repeating pattern
3062        size_t basex=0;
3063        size_t basey=0;
3064        size_t basez=0;
3065    #ifdef ESYS_MPI    
3066        basex=X*m_gNE[0]/m_NX[0];
3067        basey=Y*m_gNE[1]/m_NX[1];
3068        basez=Z*m_gNE[2]/m_NX[2];
3069            
3070    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3071        
3072    #endif    
3073        if (seed==0)
3074        {
3075        seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception
3076        }
3077        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3078    
3079    
3080    for (int z=0;z<ext[2];++z)
3081    {
3082        for (int y=0;y<ext[1];++y)
3083        {
3084        for (int x=0;x<ext[0];++x)
3085        {
3086            cout << "(";
3087            for (int p=0;p<numvals;++p)
3088            {
3089            cout << src[(x+y*ext[0]+z*ext[0]*ext[1])*numvals+p] << ", ";
3090            }
3091            cout << ") ";
3092          
3093        }
3094        cout << endl;
3095        }
3096        cout << endl;
3097    }
3098        
3099        
3100        
3101    /*
3102    cout << "Pattern:\n";    
3103    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3104    {
3105        cout << src[i++] << " ";
3106        if (i%ext[0]==0)
3107          cout << "\n";
3108        if (i%(ext[0]*ext[1])==0)
3109          cout << "\n";
3110    }*/
3111        
3112      
3113    #ifdef ESYS_MPI
3114    
3115    
3116    
3117      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);
3118      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3119            // there is an overlap of two points both of which need to have "radius" points on either side.
3120            
3121      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
3122      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3123      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3124            
3125      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);    
3126            
3127      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
3128      MPI_Status stats[50];      MPI_Status stats[50];
# Line 3019  escript::Data Brick::randomFill(long see Line 3135  escript::Data Brick::randomFill(long see
3135      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3136            
3137            
3138      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
3139    
3140        
3141    for (size_t j=0;j<incoms.size();++j)
3142    {
3143        message& m=incoms[j];
3144    for (int i=0;i<block.getBuffSize(m.srcbuffid);++i)
3145    {
3146        block.getInBuffer(m.srcbuffid)[i]=-42;
3147    }  
3148    }
3149            
3150            
3151      int comserr=0;          int comserr=0;    
# Line 3030  escript::Data Brick::randomFill(long see Line 3156  escript::Data Brick::randomFill(long see
3156      block.setUsed(m.destbuffid);      block.setUsed(m.destbuffid);
3157      }      }
3158    
3159      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3160      {      {
3161      message& m=incoms[i];      message& m=outcoms[i];
3162      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++));
3163    if (m.destID==1)
3164    {
3165    cout << "Sending to 1\n";
3166    for (int i=0;i<block.getBuffSize(m.srcbuffid);++i)
3167    {
3168        cout << block.getOutBuffer(m.srcbuffid)[i] << " ";
3169    }
3170    cout << endl;  
3171    }
3172      }          }    
3173            
3174      if (!comserr)      if (!comserr)
# Line 3050  escript::Data Brick::randomFill(long see Line 3185  escript::Data Brick::randomFill(long see
3185      }      }
3186            
3187      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3188        
3189        for (size_t i=0;i<incoms.size();++i)
3190        {
3191        message& m=incoms[i];
3192        cout << "From" <<  m.sourceID << "Recv " << (int)m.destbuffid << endl;
3193        for (int j=0;j<block.getBuffSize(m.destbuffid);++j)
3194        {
3195            cout << block.getInBuffer(m.destbuffid)[j] << " ";
3196        }
3197        cout << endl;
3198        }    
3199        
3200        
3201        
3202        
3203    
3204  #endif      #endif    
3205      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3206      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);  /*    
3207      // don't need to check for exwrite because we just made it  cout << "Pattern (post transfer):\n";    
3208      escript::DataVector& dv=resdat.getExpandedVectorReference();  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3209      double* convolution=get3DGauss(radius, sigma);  {
3210      for (size_t z=0;z<(internal[2]);++z)      cout << src[i++] << " ";
3211        if (i%ext[0]==0)
3212          cout << "\n";
3213        if (i%(ext[0]*ext[1])==0)
3214          cout << "\n";
3215    }   */
3216        
3217        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3218      {      {
3219      for (size_t y=0;y<(internal[1]);++y)            
3220        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3221        escript::Data resdat(0, shape, fs , true);
3222        // don't need to check for exwrite because we just made it
3223        escript::DataVector& dv=resdat.getExpandedVectorReference();
3224        
3225        
3226        // now we need to copy values over
3227        for (size_t z=0;z<(internal[2]);++z)
3228      {      {
3229          for (size_t x=0;x<(internal[0]);++x)          for (size_t y=0;y<(internal[1]);++y)    
3230          {              {
3231          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)
3232                    {
3233                for (unsigned int i=0;i<numvals;++i)
3234                {
3235                    dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3236                }
3237            }
3238            }
3239        }  
3240        delete[] src;
3241        return resdat;      
3242        }
3243        else        // filter enabled  
3244        {
3245        
3246        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3247        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3248        // don't need to check for exwrite because we just made it
3249        escript::DataVector& dv=resdat.getExpandedVectorReference();
3250        double* convolution=get3DGauss(radius, sigma);
3251    
3252    // cout << "Convolution matrix\n";
3253    // size_t di=(radius*2+1);
3254    // double ctot=0;
3255    // for (int i=0;i<di*di*di;++i)
3256    // {
3257    //     cout << convolution[i] << " ";
3258    //     ctot+=convolution[i];
3259    //     if ((i+1)%di==0)
3260    //     {
3261    //  cout << "\n";
3262    //     }
3263    //     if ((i+1)%(di*di)==0)
3264    //     {
3265    //  cout << "\n";
3266    //     }
3267    // }
3268    //
3269    // cout << "Sum of matrix is = " << ctot << endl;
3270    
3271        for (size_t z=0;z<(internal[2]);++z)
3272        {
3273            for (size_t y=0;y<(internal[1]);++y)    
3274            {
3275            for (size_t x=0;x<(internal[0]);++x)
3276            {    
3277                dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3278                
3279            }
3280          }          }
3281      }      }
3282        
3283    // cout << "\nResulting matrix:\n";
3284    //     for (size_t z=0;z<(internal[2]);++z)
3285    //     {
3286    //  for (size_t y=0;y<(internal[1]);++y)    
3287    //  {
3288    //      for (size_t x=0;x<(internal[0]);++x)
3289    //      {    
3290    //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";
3291    //      }
3292    //      cout << endl;
3293    //  }
3294    //  cout << endl;
3295    //     }
3296    
3297        
3298        delete[] convolution;
3299        delete[] src;
3300        return resdat;
3301        
3302      }      }
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3303  }  }
3304    
3305    
# Line 3081  escript::Data Brick::randomFill(long see Line 3310  escript::Data Brick::randomFill(long see
3310  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3311      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3312      //is the found element even owned by this rank      //is the found element even owned by this rank
3313        // (inside owned or shared elements but will map to an owned element)
3314      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3315          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3316                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3317            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3318                    + m_dx[dim]/2.;
3319            if (min > coords[dim] || max < coords[dim]) {
3320              return NOT_MINE;              return NOT_MINE;
3321          }          }
3322      }      }
# Line 3102  int Brick::findNode(const double *coords Line 3335  int Brick::findNode(const double *coords
3335          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3336      }      }
3337      //find the closest node      //find the closest node
3338      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3339          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3340          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3341              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3342              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3343                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3344                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3345                  if (total < minDist) {                  if (total < minDist) {
3346                      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],
3347                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3348                      minDist = total;                      minDist = total;
3349                  }                  }
3350              }              }
3351          }          }
3352      }      }
3353        if (closest == NOT_MINE) {
3354            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3355        }
3356      return closest;      return closest;
3357  }  }
3358    

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

  ViewVC Help
Powered by ViewVC 1.1.26