/[escript]/branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
ViewVC logotype

Contents of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3697 - (show annotations)
Mon Nov 28 04:52:00 2011 UTC (7 years, 11 months ago) by caltinay
File size: 18618 byte(s)
[x] Rectangle node id's and weipa compatibility
[ ] Brick...

1
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 pair<double,double> xdx = getFirstCoordAndSpacing(0);
142 pair<double,double> ydy = getFirstCoordAndSpacing(1);
143 #pragma omp parallel
144 {
145 #pragma omp for
146 for (dim_t i0 = 0; i0 < m_N0; i0++) {
147 coords[0][i0]=xdx.first+i0*xdx.second;
148 }
149 #pragma omp for
150 for (dim_t i1 = 0; i1 < m_N1; i1++) {
151 coords[1][i1]=ydy.first+i1*ydy.second;
152 }
153 }
154 IndexVector dims = getNumNodesPerDim();
155
156 // write mesh
157 DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
158 DB_COLLINEAR, NULL);
159
160 // 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 if (m_mpiInfo->rank == 0) {
171 vector<string> tempstrings;
172 vector<char*> names;
173 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 names.push_back((char*)tempstrings.back().c_str());
178 }
179 vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
180 DBSetDir(dbfile, "/");
181 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 }
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 const int* Rectangle::borrowSampleReferenceIDs(int fsType) const
221 {
222 switch (fsType) {
223 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 << fsType;
236 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 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 }
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 dim_t numNeighbours = 0;
326 RankVector neighbour(numNeighbours);
327 IndexVector offsetInShared(numNeighbours+1);
328 IndexVector shared(offsetInShared[numNeighbours]);
329
330 Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
331 getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
332 &offsetInShared[0], 1, 0, m_mpiInfo);
333 Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
334 getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
335 &offsetInShared[0], 1, 0, m_mpiInfo);
336 Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
337
338 // create patterns
339 dim_t M=0, N=0;
340 int* ptr=NULL;
341 index_t* index=NULL;
342 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 M, N, ptr, index);
346 Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
347 M, N, ptr, index);
348
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 return pattern;
360 */
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 pair<double,double> xdx = getFirstCoordAndSpacing(0);
372 pair<double,double> ydy = getFirstCoordAndSpacing(1);
373 for (index_t i=0; i < getNumNodes(); i++) {
374 cout << " " << setw(5) << m_nodeId[i]
375 << " " << xdx.first+(i%m_N0)*xdx.second
376 << " " << ydy.first+(i/m_N0)*ydy.second << endl;
377 }
378 }
379 }
380
381 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 //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 pair<double,double> xdx = getFirstCoordAndSpacing(0);
456 pair<double,double> ydy = getFirstCoordAndSpacing(1);
457 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 point[0] = xdx.first+i0*xdx.second;
463 point[1] = ydy.first+i1*ydy.second;
464 }
465 }
466 }
467
468 //private
469 void Rectangle::populateSampleIds()
470 {
471 // 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 m_nodeId.resize(getNumNodes());
495
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 #pragma omp parallel for
504 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 for (dim_t k=0; k<getNumNodes(); k++) {
540 cout << m_nodeId[k] << endl;
541 }
542
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 }
557
558 } // end of namespace ripley
559

  ViewVC Help
Powered by ViewVC 1.1.26