/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 4660 by sshaw, Fri Feb 7 03:08:49 2014 UTC revision 4705 by jfenwick, Fri Feb 21 02:36:15 2014 UTC
# Line 21  Line 21 
21  #include <ripley/WaveAssembler2D.h>  #include <ripley/WaveAssembler2D.h>
22  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
23  #include "esysUtils/EsysRandom.h"  #include "esysUtils/EsysRandom.h"
24    #include "blocktools.h"
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
27  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 1875  namespace Line 1876  namespace
1876      {      {
1877          double* arr=new double[(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)];
1878          double common=M_1_PI*0.5*1/(sigma*sigma);          double common=M_1_PI*0.5*1/(sigma*sigma);
1879      double total=0;          double total=0;
1880      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
1881      for (int y=-r;y<=r;++y)          for (int y=-r;y<=r;++y)
1882      {          {
1883          for (int x=-r;x<=r;++x)              for (int x=-r;x<=r;++x)
1884          {                      {        
1885              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));                  arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
1886  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;                  total+=arr[(x+r)+(y+r)*(r*2+1)];
1887              total+=arr[(x+r)+(y+r)*(r*2+1)];              }
1888          }          }
1889      }          double invtotal=1/total;
1890      double invtotal=1/total;          for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
1891  //cout << "Inv total is "    << invtotal << endl;          {
1892      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)              arr[p]*=invtotal;
1893      {          }
1894          arr[p]*=invtotal;          return arr;
 //cout << p << "->" << arr[p] << endl;        
     }  
     return arr;  
1895      }      }
1896            
1897      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 1901  namespace Line 1899  namespace
1899      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
1900      {      {
1901          size_t bx=xp-radius, by=yp-radius;          size_t bx=xp-radius, by=yp-radius;
1902      size_t sbase=bx+by*width;          size_t sbase=bx+by*width;
1903      double result=0;          double result=0;
1904      for (int y=0;y<2*radius+1;++y)          for (int y=0;y<2*radius+1;++y)
1905      {              {        
1906          for (int x=0;x<2*radius+1;++x)              for (int x=0;x<2*radius+1;++x)
1907          {              {
1908              result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];                  result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
1909          }              }
1910      }          }
1911          return result;                return result;      
1912      }      }
1913  }  }
1914    
1915    
1916    /* This is a wrapper for filtered (and non-filtered) randoms
1917     * For detailed doco see randomFillWorker
1918    */
1919    escript::Data Rectangle::randomFill(const escript::DataTypes::ShapeType& shape,
1920           const escript::FunctionSpace& what,
1921           long seed, const boost::python::tuple& filter) const
1922    {
1923        int numvals=escript::DataTypes::noValues(shape);
1924        if (len(filter)>0 && (numvals!=1))
1925        {
1926            throw RipleyException("Ripley only supports filters for scalar data.");
1927        }
1928        escript::Data res=randomFillWorker(shape, seed, filter);
1929        if (res.getFunctionSpace()!=what)
1930        {
1931            escript::Data r=escript::Data(res, what);
1932            return r;
1933        }
1934        return res;
1935    }
1936    
1937    
1938  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
1939  The dimensions of the rectangle being filled are internal[0] x internal[1] points.  The dimensions of the rectangle being filled are internal[0] x internal[1] points.
# Line 1955  that ripley has. Line 1974  that ripley has.
1974    
1975    
1976  */  */
1977  escript::Data Rectangle::randomFill(long seed, const boost::python::tuple& filter) const  escript::Data Rectangle::randomFillWorker(const escript::DataTypes::ShapeType& shape,
1978           long seed, const boost::python::tuple& filter) const
1979  {  {
1980      if (m_numDim!=2)      if (m_numDim!=2)
1981      {      {
1982          throw RipleyException("Only 2D supported at this time.");          throw RipleyException("Only 2D supported at this time.");
1983      }      }
1984      if (len(filter)!=3) {  
1985          throw RipleyException("Unsupported random filter");      unsigned int radius=0;      // these are only used by gaussian
1986      }      double sigma=0.5;
1987      boost::python::extract<string> ex(filter[0]);      
1988      if (!ex.check() || (ex()!="gaussian"))      unsigned int numvals=escript::DataTypes::noValues(shape);
1989            
1990        
1991        if (len(filter)==0)
1992      {      {
1993          throw RipleyException("Unsupported random filter");          // nothing special required here yet
1994      }      }    
1995      boost::python::extract<unsigned int> ex1(filter[1]);      else if (len(filter)==3)
     if (!ex1.check())  
1996      {      {
1997          throw RipleyException("Radius of gaussian filter must be a positive integer.");          boost::python::extract<string> ex(filter[0]);
1998            if (!ex.check() || (ex()!="gaussian"))
1999            {
2000                throw RipleyException("Unsupported random filter");
2001            }
2002            boost::python::extract<unsigned int> ex1(filter[1]);
2003            if (!ex1.check())
2004            {
2005                throw RipleyException("Radius of gaussian filter must be a positive integer.");
2006            }
2007            radius=ex1();
2008            sigma=0.5;
2009            boost::python::extract<double> ex2(filter[2]);
2010            if (!ex2.check() || (sigma=ex2())<=0)
2011            {
2012                throw RipleyException("Sigma must be a postive floating point number.");
2013            }        
2014      }      }
2015      unsigned int radius=ex1();      else
     double sigma=0.5;  
     boost::python::extract<double> ex2(filter[2]);  
     if (!ex2.check() || (sigma=ex2())<=0)  
2016      {      {
2017          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Unsupported random filter for Rectangle.");
2018      }          }
2019          
2020      
2021            
2022      size_t internal[2];      size_t internal[2];
2023      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
2024      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
2025      size_t ext[2];      size_t ext[2];
2026      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=internal[0]+2*radius;        // includes points we need as input
2027      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=internal[1]+2*radius;        // for smoothing
2028            
2029      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
2030      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
2031    
2032      if ((2*radius>=internal[0]) || (2*radius>=internal[1]))      if (2*radius>=internal[0]-4)
     {  
         throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");  
     }  
       
     
     double* src=new double[ext[0]*ext[1]];  
     esysUtils::randomFillArray(seed, src, ext[0]*ext[1]);    
       
     
 #ifdef ESYS_MPI      
     size_t inset=2*radius+1;  
     size_t Eheight=ext[1]-2*inset;  // how high the E (shared) region is  
     size_t Swidth=ext[0]-2*inset;  
   
     double* SWin=new double[inset*inset];  memset(SWin, 0, inset*inset*sizeof(double));  
     double* SEin=new double[inset*inset];  memset(SEin, 0, inset*inset*sizeof(double));  
     double* NWin=new double[inset*inset];  memset(NWin, 0, inset*inset*sizeof(double));  
     double* Sin=new double[inset*Swidth];  memset(Sin, 0, inset*Swidth*sizeof(double));  
     double* Win=new double[inset*Eheight];  memset(Win, 0, inset*Eheight*sizeof(double));  
   
     double* NEout=new double[inset*inset];  memset(NEout, 0, inset*inset*sizeof(double));  
     unsigned int base=ext[0]-inset+(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2033      {      {
2034      memcpy(NEout+inset*i, src+base, inset*sizeof(double));          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
     base+=ext[0];  
2035      }      }
2036      double* NWout=new double[inset*inset];  memset(NWout, 0, inset*inset*sizeof(double));      if (2*radius>=internal[1]-4)
     base=(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2037      {      {
2038      memcpy(NWout+inset*i, src+base, inset*sizeof(double));          throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
     base+=ext[0];  
2039      }      }
2040    
2041        double* src=new double[ext[0]*ext[1]*numvals];
2042        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*numvals);  
2043            
2044      double* SEout=new double[inset*inset];  memset(SEout, 0, inset*inset*sizeof(double));  
2045      base=ext[0]-inset;  
2046      for (int i=0;i<inset;++i)  #ifdef ESYS_MPI
2047      {      if ((internal[0]<5) || (internal[1]<5))
     memcpy(SEout+inset*i, src+base, inset*sizeof(double));  
     base+=ext[0];  
     }  
     double* Nout=new double[inset*Swidth];  memset(Nout, 0, inset*Swidth*sizeof(double));  
     base=inset+(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2048      {      {
2049      memcpy(Nout+Swidth*i, src+base, Swidth*sizeof(double));          // since the dimensions are equal for all ranks, this exception
2050      base+=ext[0];          // will be thrown on all ranks
2051            throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
2052      }      }
2053        dim_t X=m_mpiInfo->rank%m_NX[0];
2054        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2055    #endif      
2056    
2057    /*    
2058        // if we wanted to test a repeating pattern
2059        size_t basex=0;
2060        size_t basey=0;
2061    #ifdef ESYS_MPI    
2062        basex=X*m_gNE[0]/m_NX[0];
2063        basey=Y*m_gNE[1]/m_NX[1];
2064    #endif
2065            
2066        esysUtils::patternFillArray2D(ext[0], ext[1], src, 4, basex, basey, numvals);
2067    */
2068    
2069            
2070      double* Eout=new double[inset*Eheight];  memset(Eout, 0, inset*Eheight*sizeof(double));  #ifdef ESYS_MPI  
2071      base=ext[0]-inset+inset*ext[0];      
2072      for (unsigned int i=0;i<Eheight;++i)      BlockGrid2 grid(m_NX[0]-1, m_NX[1]-1);
2073      {      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
2074      memcpy(Eout+i*inset, src+base, inset*sizeof(double));                  // there is an overlap of two points both of which need to have "radius" points on either side.
2075      base+=ext[0];      
2076      }        size_t xmidlen=ext[0]-2*inset;      // how wide is the x-dimension between the two insets
2077        size_t ymidlen=ext[1]-2*inset;      
2078        
2079        Block2 block(ext[0], ext[1], inset, xmidlen, ymidlen, numvals);
2080    
2081      MPI_Request reqs[10];      MPI_Request reqs[40];               // a non-tight upper bound on how many we need
2082      MPI_Status stats[10];      MPI_Status stats[40];
2083      short rused=0;      short rused=0;
2084            
2085      dim_t X=m_mpiInfo->rank%m_NX[0];      messvec incoms;
2086      dim_t Y=m_mpiInfo->rank/m_NX[0];      messvec outcoms;  
     dim_t row=m_NX[0];  
     bool swused=false;      // These vars will be true if data needs to be copied out of them  
     bool seused=false;  
     bool nwused=false;  
     bool sused=false;  
     bool wused=false;      
2087            
2088      // Tags:      grid.generateInNeighbours(X, Y, incoms);
2089      // 10 : EW transfer (middle)      grid.generateOutNeighbours(X, Y, outcoms);
     // 8 : NS transfer (middle)  
     // 7 : NE corner -> to N, E and NE  
     // 11 : NW corner to SW corner (only used on the left hand edge  
     // 12 : SE corner to SW corner (only used on the bottom edge  
2090            
2091        block.copyAllToBuffer(src);  
   
     int comserr=0;  
     if (Y!=0)   // not on bottom row,  
     {  
     if (X!=0)   // not on the left hand edge  
     {  
         // recv bottomleft from SW  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
         comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  
         wused=true;  
     }  
     else    // on the left hand edge  
     {  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (Y-1)*row, 11, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
     }  
     comserr|=MPI_Irecv(Sin, Swidth*inset, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));  
     sused=true;  
     comserr|=MPI_Irecv(SEin, inset*inset, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     seused=true;  
   
         
     }  
     else        // on the bottom row  
     {  
     if (X!=0)  
     {  
         comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  
         wused=true;  
         // Need to use tag 12 here because SW is coming from the East not South East  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, X-1+Y*row, 12, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
     }  
     if (X!=(row-1))  
     {  
         comserr|=MPI_Isend(SEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 12, m_mpiInfo->comm, reqs+(rused++));    
     }  
     }  
2092            
2093      if (Y!=(m_NX[1]-1)) // not on the top row      int comserr=0;    
2094        for (size_t i=0;i<incoms.size();++i)
2095      {      {
2096      comserr|=MPI_Isend(Nout, inset*Swidth, MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));          message& m=incoms[i];
2097      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+(Y+1)*row, 7, 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++));
2098      if (X!=(row-1)) // not on right hand edge          block.setUsed(m.destbuffid);
2099      {      
         comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     }  
     if (X==0)   // left hand edge  
     {  
         comserr|=MPI_Isend(NWout, inset*inset, MPI_DOUBLE, (Y+1)*row,11, m_mpiInfo->comm, reqs+(rused++));        
     }    
2100      }      }
2101      if (X!=(row-1)) // not on right hand edge  
2102      {      for (size_t i=0;i<outcoms.size();++i)
     comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     comserr|=MPI_Isend(Eout, Eheight*inset, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));  
     }  
     if (X!=0)  
2103      {      {
2104      comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));          message& m=outcoms[i];
2105      nwused=true;          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
2106              }    
         
     }  
2107            
2108      if (!comserr)      if (!comserr)
2109      {          {    
2110          comserr=MPI_Waitall(rused, reqs, stats);          comserr=MPI_Waitall(rused, reqs, stats);    
2111      }      }
2112    
2113      if (comserr)      if (comserr)
2114      {      {
2115      // 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.
2116      // 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.
2117      // 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
2118          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
2119      }      }
2120            
2121        block.copyUsedFromBuffer(src);    
     // Need to copy the values back from the buffers  
     // Copy from SW  
2122            
2123      if (swused)  #endif    
2124      {      
2125      base=0;      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, SWin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (seused)  
     {  
         base=ext[0]-inset;  
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, SEin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (nwused)  
     {  
         base=(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, NWin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (sused)  
     {  
        base=inset;  
        for (unsigned int i=0;i<inset;++i)  
        {  
        memcpy(src+base, Sin+i*Swidth, Swidth*sizeof(double));  
        base+=ext[0];  
        }  
     }  
     if (wused)  
2126      {      {
     base=inset*ext[0];  
     for (unsigned int i=0;i<Eheight;++i)  
     {  
         memcpy(src+base, Win+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
2127                
2128            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2129            escript::Data resdat(0, shape, fs , true);
2130            // don't need to check for exwrite because we just made it
2131            escript::DataVector& dv=resdat.getExpandedVectorReference();
2132            
2133            
2134            // now we need to copy values over
2135            for (size_t y=0;y<(internal[1]);++y)    
2136            {
2137                for (size_t x=0;x<(internal[0]);++x)
2138                {
2139                    for (unsigned int i=0;i<numvals;++i)
2140                    {
2141                        dv[i+(x+y*(internal[0]))*numvals]=src[i+(x+y*ext[0])*numvals];
2142                    }
2143                }
2144            }
2145            delete[] src;
2146            return resdat;      
2147        }
2148        else                // filter enabled      
2149        {    
2150            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2151            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
2152            // don't need to check for exwrite because we just made it
2153            escript::DataVector& dv=resdat.getExpandedVectorReference();
2154            double* convolution=get2DGauss(radius, sigma);
2155            for (size_t y=0;y<(internal[1]);++y)    
2156            {
2157                for (size_t x=0;x<(internal[0]);++x)
2158                {    
2159                    dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);
2160                    
2161                }
2162            }
2163            delete[] convolution;
2164            delete[] src;
2165            return resdat;
2166      }      }
       
     delete[] SWin;  
     delete[] SEin;  
     delete[] NWin;  
     delete[] Sin;  
     delete[] Win;  
   
     delete[] NEout;  
     delete[] NWout;  
     delete[] SEout;  
     delete[] Nout;  
     delete[] Eout;  
 #endif      
     escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());  
     escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);  
     // don't need to check for exwrite because we just made it  
     escript::DataVector& dv=resdat.getExpandedVectorReference();  
     double* convolution=get2DGauss(radius, sigma);  
     for (size_t y=0;y<(internal[1]);++y)      
     {  
         for (size_t x=0;x<(internal[0]);++x)  
     {      
         dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);  
           
     }  
     }  
     delete[] convolution;  
     delete[] src;  
     return resdat;  
2167  }  }
2168    
2169  int Rectangle::findNode(const double *coords) const {  int Rectangle::findNode(const double *coords) const {

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

  ViewVC Help
Powered by ViewVC 1.1.26