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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3697 - (hide annotations)
Mon Nov 28 04:52:00 2011 UTC (8 years ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
File size: 18618 byte(s)
[x] Rectangle node id's and weipa compatibility
[ ] Brick...

1 caltinay 3691
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2011 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     #include <ripley/Rectangle.h>
15     extern "C" {
16     #include "paso/SystemMatrixPattern.h"
17     }
18    
19     #if USE_SILO
20     #include <silo.h>
21     #ifdef ESYS_MPI
22     #include <pmpio.h>
23     #endif
24     #endif
25    
26     #include <iomanip>
27    
28     using namespace std;
29    
30     namespace ripley {
31    
32     Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :
33     RipleyDomain(2),
34     m_gNE0(n0),
35     m_gNE1(n1),
36     m_l0(l0),
37     m_l1(l1),
38     m_NX(d0),
39     m_NY(d1)
40     {
41     // ensure number of subdivisions is valid and nodes can be distributed
42     // among number of ranks
43     if (m_NX*m_NY != m_mpiInfo->size)
44     throw RipleyException("Invalid number of spatial subdivisions");
45    
46     if (n0%m_NX > 0 || n1%m_NY > 0)
47     throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
48    
49     // local number of elements
50     m_NE0 = n0/m_NX;
51     m_NE1 = n1/m_NY;
52     // local number of nodes (not necessarily owned)
53     m_N0 = m_NE0+1;
54     m_N1 = m_NE1+1;
55     // bottom-left node is at (offset0,offset1) in global mesh
56     m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
57     m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);
58     populateSampleIds();
59     }
60    
61     Rectangle::~Rectangle()
62     {
63     }
64    
65     string Rectangle::getDescription() const
66     {
67     return "ripley::Rectangle";
68     }
69    
70     bool Rectangle::operator==(const AbstractDomain& other) const
71     {
72     if (dynamic_cast<const Rectangle*>(&other))
73     return this==&other;
74    
75     return false;
76     }
77    
78     void Rectangle::dump(const string& fileName) const
79     {
80     #if USE_SILO
81     string fn(fileName);
82     if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
83     fn+=".silo";
84     }
85    
86     const int NUM_SILO_FILES = 1;
87     const char* blockDirFmt = "/block%04d";
88     int driver=DB_HDF5;
89     string siloPath;
90     DBfile* dbfile = NULL;
91    
92     #ifdef ESYS_MPI
93     PMPIO_baton_t* baton = NULL;
94     #endif
95    
96     if (m_mpiInfo->size > 1) {
97     #ifdef ESYS_MPI
98     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
99     0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
100     PMPIO_DefaultClose, (void*)&driver);
101     // try the fallback driver in case of error
102     if (!baton && driver != DB_PDB) {
103     driver = DB_PDB;
104     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
105     0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
106     PMPIO_DefaultClose, (void*)&driver);
107     }
108     if (baton) {
109     char str[64];
110     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
111     siloPath = str;
112     dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
113     }
114     #endif
115     } else {
116     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
117     getDescription().c_str(), driver);
118     // try the fallback driver in case of error
119     if (!dbfile && driver != DB_PDB) {
120     driver = DB_PDB;
121     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
122     getDescription().c_str(), driver);
123     }
124     }
125    
126     if (!dbfile)
127     throw RipleyException("dump: Could not create Silo file");
128    
129     /*
130     if (driver==DB_HDF5) {
131     // gzip level 1 already provides good compression with minimal
132     // performance penalty. Some tests showed that gzip levels >3 performed
133     // rather badly on escript data both in terms of time and space
134     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
135     }
136     */
137    
138     boost::scoped_ptr<double> x(new double[m_N0]);
139     boost::scoped_ptr<double> y(new double[m_N1]);
140     double* coords[2] = { x.get(), y.get() };
141 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
142     pair<double,double> ydy = getFirstCoordAndSpacing(1);
143 caltinay 3691 #pragma omp parallel
144     {
145     #pragma omp for
146     for (dim_t i0 = 0; i0 < m_N0; i0++) {
147 caltinay 3697 coords[0][i0]=xdx.first+i0*xdx.second;
148 caltinay 3691 }
149     #pragma omp for
150     for (dim_t i1 = 0; i1 < m_N1; i1++) {
151 caltinay 3697 coords[1][i1]=ydy.first+i1*ydy.second;
152 caltinay 3691 }
153     }
154 caltinay 3697 IndexVector dims = getNumNodesPerDim();
155    
156     // write mesh
157     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
158 caltinay 3691 DB_COLLINEAR, NULL);
159    
160 caltinay 3697 // write node ids
161     DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,
162     NULL, 0, DB_INT, DB_NODECENT, NULL);
163    
164     // write element ids
165     dims = getNumElementsPerDim();
166     DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
167     &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
168    
169     // rank 0 writes multimesh and multivar
170 caltinay 3691 if (m_mpiInfo->rank == 0) {
171     vector<string> tempstrings;
172 caltinay 3697 vector<char*> names;
173 caltinay 3691 for (dim_t i=0; i<m_mpiInfo->size; i++) {
174     stringstream path;
175     path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
176     tempstrings.push_back(path.str());
177 caltinay 3697 names.push_back((char*)tempstrings.back().c_str());
178 caltinay 3691 }
179 caltinay 3697 vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
180 caltinay 3691 DBSetDir(dbfile, "/");
181 caltinay 3697 DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
182     &types[0], NULL);
183     tempstrings.clear();
184     names.clear();
185     for (dim_t i=0; i<m_mpiInfo->size; i++) {
186     stringstream path;
187     path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
188     tempstrings.push_back(path.str());
189     names.push_back((char*)tempstrings.back().c_str());
190     }
191     types.assign(m_mpiInfo->size, DB_QUADVAR);
192     DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
193     &types[0], NULL);
194     tempstrings.clear();
195     names.clear();
196     for (dim_t i=0; i<m_mpiInfo->size; i++) {
197     stringstream path;
198     path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
199     tempstrings.push_back(path.str());
200     names.push_back((char*)tempstrings.back().c_str());
201     }
202     DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
203     &types[0], NULL);
204 caltinay 3691 }
205    
206     if (m_mpiInfo->size > 1) {
207     #ifdef ESYS_MPI
208     PMPIO_HandOffBaton(baton, dbfile);
209     PMPIO_Finish(baton);
210     #endif
211     } else {
212     DBClose(dbfile);
213     }
214    
215     #else // USE_SILO
216     throw RipleyException("dump(): no Silo support");
217     #endif
218     }
219    
220 caltinay 3697 const int* Rectangle::borrowSampleReferenceIDs(int fsType) const
221 caltinay 3691 {
222 caltinay 3697 switch (fsType) {
223 caltinay 3691 case Nodes:
224     return &m_nodeId[0];
225     case Elements:
226     return &m_elementId[0];
227     case FaceElements:
228     return &m_faceId[0];
229     default:
230     break;
231     }
232    
233     stringstream msg;
234     msg << "borrowSampleReferenceIDs() not implemented for function space type "
235 caltinay 3697 << fsType;
236 caltinay 3691 throw RipleyException(msg.str());
237     }
238    
239     bool Rectangle::ownSample(int fsCode, index_t id) const
240     {
241     #ifdef ESYS_MPI
242     if (fsCode == Nodes) {
243     const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
244     const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
245     return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
246     } else
247     throw RipleyException("ownSample() only implemented for Nodes");
248     #else
249     return true;
250     #endif
251     }
252    
253     void Rectangle::interpolateOnDomain(escript::Data& target,
254     const escript::Data& in) const
255     {
256     const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));
257     const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));
258     if (inDomain != *this)
259     throw RipleyException("Illegal domain of interpolant");
260     if (targetDomain != *this)
261     throw RipleyException("Illegal domain of interpolation target");
262    
263 caltinay 3697 stringstream msg;
264     msg << "interpolateOnDomain() not implemented for function space "
265     << in.getFunctionSpace().getTypeCode() << " -> "
266     << target.getFunctionSpace().getTypeCode();
267    
268     switch (in.getFunctionSpace().getTypeCode()) {
269     case Nodes:
270     case DegreesOfFreedom:
271     switch (target.getFunctionSpace().getTypeCode()) {
272     case Elements:
273     {
274     const double tmp0_2 = 0.62200846792814621559;
275     const double tmp0_1 = 0.044658198738520451079;
276     const double tmp0_0 = 0.16666666666666666667;
277     const dim_t numComp = in.getDataPointSize();
278     escript::Data* inn=const_cast<escript::Data*>(&in);
279     #pragma omp parallel for
280     for (index_t k1=0; k1 < m_NE1; ++k1) {
281     for (index_t k0=0; k0 < m_NE0; ++k0) {
282     const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));
283     const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
284     const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));
285     const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));
286     double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));
287     for (index_t i=0; i < numComp; ++i) {
288     o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
289     o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
290     o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
291     o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
292     } // close component loop i
293     } // close k0 loop
294     } // close k1 loop
295     break;
296     }
297    
298     case DegreesOfFreedom:
299     target=in;
300     break;
301    
302     default:
303     throw RipleyException(msg.str());
304     }
305     default:
306     throw RipleyException(msg.str());
307     }
308 caltinay 3691 }
309    
310     Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
311     bool reducedColOrder) const
312     {
313     if (reducedRowOrder || reducedColOrder)
314     throw RipleyException("getPattern() not implemented for reduced order");
315    
316     /*
317     // create distribution
318     IndexVector dist;
319     for (index_t i=0; i<m_mpiInfo->size+1; i++)
320     dist.push_back(i*getNumNodes());
321     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
322     &dist[0], 1, 0);
323    
324     // connectors
325 caltinay 3697 dim_t numNeighbours = 0;
326 caltinay 3691 RankVector neighbour(numNeighbours);
327     IndexVector offsetInShared(numNeighbours+1);
328     IndexVector shared(offsetInShared[numNeighbours]);
329    
330     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
331 caltinay 3697 getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
332     &offsetInShared[0], 1, 0, m_mpiInfo);
333 caltinay 3691 Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
334 caltinay 3697 getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
335     &offsetInShared[0], 1, 0, m_mpiInfo);
336 caltinay 3691 Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
337    
338     // create patterns
339 caltinay 3697 dim_t M=0, N=0;
340     int* ptr=NULL;
341     index_t* index=NULL;
342 caltinay 3691 Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
343     M, N, ptr, index);
344     Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
345 caltinay 3697 M, N, ptr, index);
346 caltinay 3691 Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
347 caltinay 3697 M, N, ptr, index);
348 caltinay 3691
349     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
350     MATRIX_FORMAT_DEFAULT, distribution, distribution,
351     mainPattern, colCouplePattern, rowCouplePattern,
352     connector, connector);
353     Paso_Pattern_free(mainPattern);
354     Paso_Pattern_free(colCouplePattern);
355     Paso_Pattern_free(rowCouplePattern);
356     Paso_Distribution_free(distribution);
357     Paso_SharedComponents_free(snd_shcomp);
358     Paso_SharedComponents_free(rcv_shcomp);
359 caltinay 3697 return pattern;
360 caltinay 3691 */
361     throw RipleyException("getPattern() not implemented");
362     }
363    
364     void Rectangle::Print_Mesh_Info(const bool full) const
365     {
366     RipleyDomain::Print_Mesh_Info(full);
367     if (full) {
368     cout << " Id Coordinates" << endl;
369     cout.precision(15);
370     cout.setf(ios::scientific, ios::floatfield);
371 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
372     pair<double,double> ydy = getFirstCoordAndSpacing(1);
373 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
374     cout << " " << setw(5) << m_nodeId[i]
375 caltinay 3697 << " " << xdx.first+(i%m_N0)*xdx.second
376     << " " << ydy.first+(i/m_N0)*ydy.second << endl;
377 caltinay 3691 }
378     }
379     }
380    
381 caltinay 3697 IndexVector Rectangle::getNumNodesPerDim() const
382     {
383     IndexVector ret;
384     ret.push_back(m_N0);
385     ret.push_back(m_N1);
386     return ret;
387     }
388    
389     IndexVector Rectangle::getNumElementsPerDim() const
390     {
391     IndexVector ret;
392     ret.push_back(m_NE0);
393     ret.push_back(m_NE1);
394     return ret;
395     }
396    
397     IndexVector Rectangle::getNumFacesPerBoundary() const
398     {
399     IndexVector ret(4, 0);
400     //left
401     if (m_offset0==0)
402     ret[0]=m_NE1;
403     //right
404     if (m_mpiInfo->rank%m_NX==m_NX-1)
405     ret[1]=m_NE1;
406     //bottom
407     if (m_offset1==0)
408     ret[2]=m_NE0;
409     //top
410     if (m_mpiInfo->rank/m_NX==m_NY-1)
411     ret[3]=m_NE0;
412     return ret;
413     }
414    
415     pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
416     {
417     if (dim==0) {
418     return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
419     } else if (dim==1) {
420     return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
421     }
422     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
423     }
424    
425 caltinay 3691 //protected
426     dim_t Rectangle::getNumFaceElements() const
427     {
428     dim_t n=0;
429     //left
430     if (m_offset0==0)
431     n+=m_NE1;
432     //right
433     if (m_mpiInfo->rank%m_NX==m_NX-1)
434     n+=m_NE1;
435     //bottom
436     if (m_offset1==0)
437     n+=m_NE0;
438     //top
439     if (m_mpiInfo->rank/m_NX==m_NY-1)
440     n+=m_NE0;
441    
442     return n;
443     }
444    
445     //protected
446     void Rectangle::assembleCoordinates(escript::Data& arg) const
447     {
448     escriptDataC x = arg.getDataC();
449     int numDim = m_numDim;
450     if (!isDataPointShapeEqual(&x, 1, &numDim))
451     throw RipleyException("setToX: Invalid Data object shape");
452     if (!numSamplesEqual(&x, 1, getNumNodes()))
453     throw RipleyException("setToX: Illegal number of samples in Data object");
454    
455 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
456     pair<double,double> ydy = getFirstCoordAndSpacing(1);
457 caltinay 3691 arg.requireWrite();
458     #pragma omp parallel for
459     for (dim_t i1 = 0; i1 < m_N1; i1++) {
460     for (dim_t i0 = 0; i0 < m_N0; i0++) {
461     double* point = arg.getSampleDataRW(i0+m_N0*i1);
462 caltinay 3697 point[0] = xdx.first+i0*xdx.second;
463     point[1] = ydy.first+i1*ydy.second;
464 caltinay 3691 }
465     }
466     }
467    
468     //private
469     void Rectangle::populateSampleIds()
470     {
471 caltinay 3697 // identifiers are ordered from left to right, bottom to top on each rank,
472     // except for the shared nodes which are owned by the rank below / to the
473     // left of the current rank
474    
475     // build node distribution vector first.
476     // m_nodeDistribution[i] is the first node id on rank i, that is
477     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
478     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
479     m_nodeDistribution[1]=getNumNodes();
480     for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
481     const index_t x=k%m_NX;
482     const index_t y=k/m_NX;
483     index_t numNodes=getNumNodes();
484     if (x>0)
485     numNodes-=m_N1;
486     if (y>0)
487     numNodes-=m_N0;
488     if (x>0 && y>0)
489     numNodes++; // subtracted corner twice -> fix that
490     m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
491     }
492     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
493    
494 caltinay 3691 m_nodeId.resize(getNumNodes());
495 caltinay 3697
496     // the bottom row and left column are not owned by this rank so the
497     // identifiers need to be computed accordingly
498     const index_t left = (m_offset0==0 ? 0 : 1);
499     const index_t bottom = (m_offset1==0 ? 0 : 1);
500     if (left>0) {
501     const int neighbour=m_mpiInfo->rank-1;
502     const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
503 caltinay 3691 #pragma omp parallel for
504 caltinay 3697 for (dim_t i1=bottom; i1<m_N1; i1++) {
505     m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
506     + (i1-bottom+1)*leftN0-1;
507     }
508     }
509     if (bottom>0) {
510     const int neighbour=m_mpiInfo->rank-m_NX;
511     const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
512     const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
513     #pragma omp parallel for
514     for (dim_t i0=left; i0<m_N0; i0++) {
515     m_nodeId[i0]=m_nodeDistribution[neighbour]
516     + (bottomN1-1)*bottomN0 + i0 - left;
517     }
518     }
519     if (left>0 && bottom>0) {
520     const int neighbour=m_mpiInfo->rank-m_NX-1;
521     const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
522     const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
523     m_nodeId[0]=m_nodeDistribution[neighbour]+bottomN1*bottomN0-1;
524     }
525    
526     // the rest of the id's are contiguous
527     const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
528     #pragma omp parallel for
529     for (dim_t i1=bottom; i1<m_N1; i1++) {
530     for (dim_t i0=left; i0<m_N0; i0++) {
531     m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);
532     }
533     }
534    
535     for (dim_t k=0; k<m_nodeDistribution.size(); k++) {
536     cout << m_nodeDistribution[k] << endl;
537     }
538     cout<<endl;
539 caltinay 3691 for (dim_t k=0; k<getNumNodes(); k++) {
540 caltinay 3697 cout << m_nodeId[k] << endl;
541 caltinay 3691 }
542 caltinay 3697
543     // elements
544     m_elementId.resize(getNumElements());
545     #pragma omp parallel for
546     for (dim_t k=0; k<getNumElements(); k++) {
547     m_elementId[k]=k;
548     }
549    
550     // face elements
551     m_faceId.resize(getNumFaceElements());
552     #pragma omp parallel for
553     for (dim_t k=0; k<getNumFaceElements(); k++) {
554     m_faceId[k]=k;
555     }
556 caltinay 3691 }
557    
558     } // end of namespace ripley
559    

  ViewVC Help
Powered by ViewVC 1.1.26