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

revision 4703 by jfenwick, Fri Feb 21 00:55:31 2014 UTC revision 4765 by sshaw, Wed Mar 19 00:17:16 2014 UTC
# Line 14  Line 14
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16
17    #include <algorithm>
18
19  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
20  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
21  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
22  #include <ripley/DefaultAssembler2D.h>  #include <ripley/DefaultAssembler2D.h>
23  #include <ripley/WaveAssembler2D.h>  #include <ripley/WaveAssembler2D.h>
24    #include <ripley/LameAssembler2D.h>
25    #include <ripley/domainhelpers.h>
26  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
27  #include "esysUtils/EsysRandom.h"  #include "esysUtils/EsysRandom.h"
28  #include "blocktools.h"  #include "blocktools.h"
# Line 55  Rectangle::Rectangle(int n0, int n1, dou Line 59  Rectangle::Rectangle(int n0, int n1, dou
59      }      }
60
61      bool warn=false;      bool warn=false;
62      // if number of subdivisions is non-positive, try to subdivide by the same      std::vector<int> factors;
63      // ratio as the number of elements      int ranks = m_mpiInfo->size;
64      if (d0<=0 && d1<=0) {      int epr[2] = {n0,n1};
65          warn=true;      int d[2] = {d0,d1};
66          d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));      if (d0<=0 || d1<=0) {
67          d1=m_mpiInfo->size/d0;          for (int i = 0; i < 2; i++) {
68          if (d0*d1 != m_mpiInfo->size) {              if (d[i] < 1) {
69              // ratios not the same so subdivide side with more elements only                  d[i] = 1;
70              if (n0>n1) {                  continue;
d0=0;
d1=1;
} else {
d0=1;
d1=0;
71              }              }
72                epr[i] = -1; // can no longer be max
73                //remove
74                if (ranks % d[i] != 0) {
75                    throw RipleyException("Invalid number of spatial subdivisions");
76                }
77                ranks /= d[i];
78            }
79            factorise(factors, ranks);
80            if (factors.size() != 0) {
81                warn = true;
82          }          }
83      }      }
84      if (d0<=0) {
85          // d1 is preset, determine d0 - throw further down if result is no good      while (factors.size() > 0) {
86          d0=m_mpiInfo->size/d1;          int i = epr[0] > epr[1] ? 0 : 1;
87      } else if (d1<=0) {          int f = factors.back();
88          // d0 is preset, determine d1 - throw further down if result is no good          factors.pop_back();
89          d1=m_mpiInfo->size/d0;          d[i] *= f;
90            epr[i] /= f;
91      }      }
92        d0 = d[0]; d1 = d[1];
93
94      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
95      // among number of ranks      // among number of ranks
96      if (d0*d1 != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
333      }      }
334  }  }
335
336    void Rectangle::readBinaryGridFromZipped(escript::Data& out, string filename,
338    {
339        // the mapping is not universally correct but should work on our
340        // supported platforms
341        switch (params.dataType) {
342            case DATATYPE_INT32:
344                break;
345            case DATATYPE_FLOAT32:
347                break;
348            case DATATYPE_FLOAT64:
350                break;
351            default:
352                throw RipleyException("readBinaryGridFromZipped(): invalid or unsupported datatype");
353        }
354    }
355
356  template<typename ValueType>  template<typename ValueType>
357  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
441      f.close();      f.close();
442  }  }
443
444    template<typename ValueType>
445    void Rectangle::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
447    {
448        // check destination function space
449        int myN0, myN1;
450        if (out.getFunctionSpace().getTypeCode() == Nodes) {
451            myN0 = m_NN[0];
452            myN1 = m_NN[1];
453        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
454                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
455            myN0 = m_NE[0];
456            myN1 = m_NE[1];
457        } else
458            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
459
460        // check file existence and size
461        ifstream f(filename.c_str(), ifstream::binary);
462        if (f.fail()) {
463            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
464        }
465        f.seekg(0, ios::end);
466        const int numComp = out.getDataPointSize();
467        int filesize = f.tellg();
468        f.seekg(0, ios::beg);
469        std::vector<char> compressed(filesize);
471        f.close();
472        std::vector<char> decompressed = unzip(compressed);
473        filesize = decompressed.size();
474        const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
475        if (filesize < reqsize) {
476            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
477        }
478
479        // check if this rank contributes anything
480        if (params.first[0] >= m_offset[0]+myN0 ||
481                params.first[0]+params.numValues[0] <= m_offset[0] ||
482                params.first[1] >= m_offset[1]+myN1 ||
483                params.first[1]+params.numValues[1] <= m_offset[1]) {
484            f.close();
485            return;
486        }
487
488        // now determine how much this rank has to write
489
490        // first coordinates in data object to write to
491        const int first0 = max(0, params.first[0]-m_offset[0]);
492        const int first1 = max(0, params.first[1]-m_offset[1]);
493        // indices to first value in file
494        const int idx0 = max(0, m_offset[0]-params.first[0]);
495        const int idx1 = max(0, m_offset[1]-params.first[1]);
496        // number of values to read
497        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
498        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
499
500        out.requireWrite();
501        vector<ValueType> values(num0*numComp);
502        const int dpp = out.getNumDataPointsPerSample();
503
504        for (int y=0; y<num1; y++) {
505            const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
506                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
507            for (int x=0; x<num0; x++) {
508                const int baseIndex = first0+x*params.multiplier[0]
509                                        +(first1+y*params.multiplier[1])*myN0;
510                for (int m1=0; m1<params.multiplier[1]; m1++) {
511                    for (int m0=0; m0<params.multiplier[0]; m0++) {
512                        const int dataIndex = baseIndex+m0+m1*myN0;
513                        double* dest = out.getSampleDataRW(dataIndex);
514                        for (int c=0; c<numComp; c++) {
515                            ValueType val = values[x*numComp+c];
516
517                            if (params.byteOrder != BYTEORDER_NATIVE) {
518                                char* cval = reinterpret_cast<char*>(&val);
519                                // this will alter val!!
520                                byte_swap32(cval);
521                            }
522                            if (!std::isnan(val)) {
523                                for (int q=0; q<dpp; q++) {
524                                    *dest++ = static_cast<double>(val);
525                                }
526                            }
527                        }
528                    }
529                }
530            }
531        }
532
533        f.close();
534    }
535
536  void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,  void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
537                                  int byteOrder, int dataType) const                                  int byteOrder, int dataType) const
538  {  {
# Line 1512  void Rectangle::createPattern() Line 1635  void Rectangle::createPattern()
1635      RankVector neighbour;      RankVector neighbour;
1636      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1637      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1638      int numShared=0;      int numShared=0, expectedShared=0;
1639      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
1640      const int y=m_mpiInfo->rank/m_NX[0];      const int y=m_mpiInfo->rank/m_NX[0];
1641        if (x > 0)
1642            expectedShared += nDOF1;
1643        if (x < m_NX[0] - 1)
1644            expectedShared += nDOF1;
1645        if (y > 0)
1646            expectedShared += nDOF0;
1647        if (y < m_NX[1] - 1)
1648            expectedShared += nDOF0;
1649        if (x > 0 && y > 0) expectedShared++;
1650        if (x > 0 && y < m_NX[1] - 1) expectedShared++;
1651        if (x < m_NX[0] - 1 && y > 0) expectedShared++;
1652        if (x < m_NX[0] - 1 && y < m_NX[1] - 1) expectedShared++;
1653
1654        vector<IndexVector> rowIndices(expectedShared);
1655
1656      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1657          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1658              // skip this rank              // skip this rank
# Line 1534  void Rectangle::createPattern() Line 1672  void Rectangle::createPattern()
1672                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
1673                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
1674                          if (i>0)                          if (i>0)
1675                              colIndices[firstDOF+i-1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
1676                          colIndices[firstDOF+i].push_back(numShared);                          doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
1677                          if (i<nDOF0-1)                          if (i<nDOF0-1)
1678                              colIndices[firstDOF+i+1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
1679                          m_dofMap[firstNode+i]=numDOF+numShared;                          m_dofMap[firstNode+i]=numDOF+numShared;
1680                      }                      }
1681                  } else if (i1==0) {                  } else if (i1==0) {
# Line 1549  void Rectangle::createPattern() Line 1687  void Rectangle::createPattern()
1687                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
1688                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
1689                          if (i>0)                          if (i>0)
1690                              colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
1691                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
1692                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1693                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
1694                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1695                      }                      }
1696                  } else {                  } else {
# Line 1562  void Rectangle::createPattern() Line 1700  void Rectangle::createPattern()
1700                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1701                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1702                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
1703                      colIndices[dof].push_back(numShared);                      doublyLink(colIndices, rowIndices, dof, numShared);
1704                      m_dofMap[node]=numDOF+numShared;                      m_dofMap[node]=numDOF+numShared;
1705                      ++numShared;                      ++numShared;
1706                  }                  }
1707              }              }
1708          }          }
1709      }      }
1710
1711    #pragma omp parallel for
1712        for (int i = 0; i < numShared; i++) {
1713            std::sort(rowIndices[i].begin(), rowIndices[i].end());
1714        }
1715
1716      // create connector      // create connector
1717      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1718              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 1584  void Rectangle::createPattern() Line 1727  void Rectangle::createPattern()
1727      // create main and couple blocks      // create main and couple blocks
1728      Paso_Pattern *mainPattern = createMainPattern();      Paso_Pattern *mainPattern = createMainPattern();
1729      Paso_Pattern *colPattern, *rowPattern;      Paso_Pattern *colPattern, *rowPattern;
1730      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
1731
1732      // allocate paso distribution      // allocate paso distribution
1733      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
# Line 1883  namespace Line 2026  namespace
2026              for (int x=-r;x<=r;++x)              for (int x=-r;x<=r;++x)
2027              {                      {
2028                  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));
// cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;
cerr << "arr[" << (x+r)+(y+r)*(r*2+1) << "] = " << arr[(x+r)+(y+r)*(r*2+1)]<< " ";
2029                  total+=arr[(x+r)+(y+r)*(r*2+1)];                  total+=arr[(x+r)+(y+r)*(r*2+1)];
2030              }              }
cerr << endl;
2031          }          }
cerr << "Total is " << total << endl;
2032          double invtotal=1/total;          double invtotal=1/total;
cerr << "Inv total is "        << invtotal << endl;
2034          {          {
2035              arr[p]*=invtotal;              arr[p]*=invtotal;
//cout << p << "->" << arr[p] << endl;
2036          }          }
2037          return arr;          return arr;
2038      }      }
# Line 2060  escript::Data Rectangle::randomFillWorke Line 2197  escript::Data Rectangle::randomFillWorke
2197      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];
2198  #endif        #endif
2199
2200        /*
2201      // if we wanted to test a repeating pattern      // if we wanted to test a repeating pattern
2202      size_t basex=0;      size_t basex=0;
2203      size_t basey=0;      size_t basey=0;
# Line 2070  escript::Data Rectangle::randomFillWorke Line 2207  escript::Data Rectangle::randomFillWorke
2207  #endif  #endif
2208
2209      esysUtils::patternFillArray2D(ext[0], ext[1], src, 4, basex, basey, numvals);      esysUtils::patternFillArray2D(ext[0], ext[1], src, 4, basex, basey, numvals);
2210    */
2211
cout << "Initial pattern\n";
for (int p=0;p<ext[1];++p)
{
for (int q=0;q<ext[0];++q)
{
cout << '(';
for (int d=0;d<numvals;++d)
{
cout << src[(q+p*ext[0])*numvals+d] << " ";
}
cout << ") ";
}
cout << endl;
}

2212
2213  #ifdef ESYS_MPI    #ifdef ESYS_MPI
2214
# Line 2095  for (int p=0;p<ext[1];++p) Line 2218  for (int p=0;p<ext[1];++p)
2218
2219      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
2220      size_t ymidlen=ext[1]-2*inset;            size_t ymidlen=ext[1]-2*inset;

cerr << "xmidlen=" << xmidlen << " ymidlen=" << ymidlen << endl;

2221
2222      Block2 block(ext[0], ext[1], inset, xmidlen, ymidlen, numvals);      Block2 block(ext[0], ext[1], inset, xmidlen, ymidlen, numvals);
2223
# Line 2113  cerr << "xmidlen=" << xmidlen << " ymidl Line 2232  cerr << "xmidlen=" << xmidlen << " ymidl
2232      grid.generateOutNeighbours(X, Y, outcoms);      grid.generateOutNeighbours(X, Y, outcoms);
2233
2234      block.copyAllToBuffer(src);        block.copyAllToBuffer(src);

// for (int i=0;i<9;++i)
// {
//     if (i!=4)
//     {
//       for (int j=0;j<block.getBuffSize(i);++j)
//       {
//    block.getOutBuffer(i)[j]=100+i;
//       }
//     }
// }

for (size_t j=0;j<incoms.size();++j)
{
message& m=incoms[j];
for (int i=0;i<block.getBuffSize(m.srcbuffid);++i)
{
block.getInBuffer(m.srcbuffid)[i]=-42;
}
}

2235
2236      int comserr=0;          int comserr=0;
2237      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
# Line 2144  for (int i=0;i<block.getBuffSize(m.srcbu Line 2239  for (int i=0;i<block.getBuffSize(m.srcbu
2239          message& m=incoms[i];          message& m=incoms[i];
2240          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++));
2241          block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);

2242      }      }
2243
2244      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
2245      {      {
2246          message& m=outcoms[i];          message& m=outcoms[i];
cout << "Sending " <<  (int)m.srcbuffid << " with tag " << m.tag << endl;
2247          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++));
for (int i=0;i<block.getBuffSize(m.srcbuffid);++i)
{
cout << block.getOutBuffer(m.srcbuffid)[i] << " ";
}
cout << endl;
2248      }          }
2249
2250      if (!comserr)      if (!comserr)
2251      {          {
2252          comserr=MPI_Waitall(rused, reqs, stats);          comserr=MPI_Waitall(rused, reqs, stats);

for (size_t i=0;i<incoms.size();++i)
{
message& m=incoms[i];
cout << "Gettinging " <<    (int)m.destbuffid << " with tag " << m.tag << endl;
for (int i=0;i<block.getBuffSize(m.destbuffid);++i)
{
cout << block.getInBuffer(m.destbuffid)[i] << " ";
}
cout << endl;

}

2253      }      }
2254
2255      if (comserr)      if (comserr)
# Line 2185  cout << endl; Line 2259  cout << endl;
2259          // 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
2260          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");
2261      }      }
/*
for (int i=0;i<9;++i)
{
if (i!=4)
{
for (int j=0;j<block.getBuffSize(i);++j)
{
block.getInBuffer(i)[j]=200+i;
}
}
}  */

2262
2263      block.copyUsedFromBuffer(src);          block.copyUsedFromBuffer(src);
2264
# Line 2294  int Rectangle::findNode(const double *co Line 2356  int Rectangle::findNode(const double *co
2356  void Rectangle::setAssembler(std::string type, std::map<std::string,  void Rectangle::setAssembler(std::string type, std::map<std::string,
2357          escript::Data> constants) {          escript::Data> constants) {
2358      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
2359            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
2360                throw RipleyException("Domain already using a different custom assembler");
2361            assembler_type = WAVE_ASSEMBLER;
2362          delete assembler;          delete assembler;
2363          assembler = new WaveAssembler2D(this, m_dx, m_NX, m_NE, m_NN, constants);          assembler = new WaveAssembler2D(this, m_dx, m_NX, m_NE, m_NN, constants);
2364        } else if (type.compare("LameAssembler") == 0) {
2365            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
2366                throw RipleyException("Domain already using a different custom assembler");
2367            assembler_type = LAME_ASSEMBLER;
2368            delete assembler;
2369            assembler = new LameAssembler2D(this, m_dx, m_NX, m_NE, m_NN);
2370      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
2371          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Rectangle does not support the"
2372                                  " requested assembler");                                  " requested assembler");

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