/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Annotation of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3697 - (hide annotations)
Mon Nov 28 04:52:00 2011 UTC (7 years, 3 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 11221 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/Brick.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     Brick::Brick(int n0, int n1, int n2, double l0, double l1, double l2, int d0,
33     int d1, int d2) :
34     RipleyDomain(3),
35     m_gNE0(n0),
36     m_gNE1(n1),
37     m_gNE2(n2),
38     m_l0(l0),
39     m_l1(l1),
40     m_l2(l2),
41     m_NX(d0),
42     m_NY(d1),
43     m_NZ(d2)
44     {
45     // ensure number of subdivisions is valid and nodes can be distributed
46     // among number of ranks
47     if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
48     throw RipleyException("Invalid number of spatial subdivisions");
49    
50     if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)
51     throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
52    
53     // local number of elements
54     m_NE0 = n0/m_NX;
55     m_NE1 = n1/m_NY;
56     m_NE2 = n2/m_NZ;
57     // local number of nodes (not necessarily owned)
58     m_N0 = m_NE0+1;
59     m_N1 = m_NE1+1;
60     m_N2 = m_NE2+1;
61     // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
62     m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
63     m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);
64     m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));
65     populateSampleIds();
66     }
67    
68    
69     Brick::~Brick()
70     {
71     }
72    
73     string Brick::getDescription() const
74     {
75     return "ripley::Brick";
76     }
77    
78     bool Brick::operator==(const AbstractDomain& other) const
79     {
80     if (dynamic_cast<const Brick*>(&other))
81     return this==&other;
82    
83     return false;
84     }
85    
86     void Brick::dump(const string& fileName) const
87     {
88     #if USE_SILO
89     string fn(fileName);
90     if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
91     fn+=".silo";
92     }
93    
94     const int NUM_SILO_FILES = 1;
95     const char* blockDirFmt = "/block%04d";
96     int driver=DB_HDF5;
97     string siloPath;
98     DBfile* dbfile = NULL;
99    
100     #ifdef ESYS_MPI
101     PMPIO_baton_t* baton = NULL;
102     #endif
103    
104     if (m_mpiInfo->size > 1) {
105     #ifdef ESYS_MPI
106     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
107     0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
108     PMPIO_DefaultClose, (void*)&driver);
109     // try the fallback driver in case of error
110     if (!baton && driver != DB_PDB) {
111     driver = DB_PDB;
112     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
113     0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
114     PMPIO_DefaultClose, (void*)&driver);
115     }
116     if (baton) {
117     char str[64];
118     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
119     siloPath = str;
120     dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
121     }
122     #endif
123     } else {
124     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
125     getDescription().c_str(), driver);
126     // try the fallback driver in case of error
127     if (!dbfile && driver != DB_PDB) {
128     driver = DB_PDB;
129     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
130     getDescription().c_str(), driver);
131     }
132     }
133    
134     if (!dbfile)
135     throw RipleyException("dump: Could not create Silo file");
136    
137     /*
138     if (driver==DB_HDF5) {
139     // gzip level 1 already provides good compression with minimal
140     // performance penalty. Some tests showed that gzip levels >3 performed
141     // rather badly on escript data both in terms of time and space
142     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
143     }
144     */
145    
146     boost::scoped_ptr<double> x(new double[m_N0]);
147     boost::scoped_ptr<double> y(new double[m_N1]);
148     boost::scoped_ptr<double> z(new double[m_N2]);
149     double* coords[3] = { x.get(), y.get(), z.get() };
150     #pragma omp parallel
151     {
152     #pragma omp for
153     for (dim_t i0 = 0; i0 < m_N0; i0++) {
154     coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;
155     }
156     #pragma omp for
157     for (dim_t i1 = 0; i1 < m_N1; i1++) {
158     coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;
159     }
160     #pragma omp for
161     for (dim_t i2 = 0; i2 < m_N2; i2++) {
162     coords[2][i2]=(m_l2*(i2+m_offset2))/m_gNE2;
163     }
164     }
165     int dims[] = { m_N0, m_N1, m_N2 };
166     DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,
167     DB_COLLINEAR, NULL);
168    
169 caltinay 3697 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 3, NULL,
170     0, DB_INT, DB_NODECENT, NULL);
171 caltinay 3691
172     if (m_mpiInfo->rank == 0) {
173     vector<string> tempstrings;
174     vector<char*> names;
175     for (dim_t i=0; i<m_mpiInfo->size; i++) {
176     stringstream path;
177     path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
178     tempstrings.push_back(path.str());
179     names.push_back((char*)tempstrings.back().c_str());
180     }
181     vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
182     DBSetDir(dbfile, "/");
183     DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
184     &types[0], NULL);
185     tempstrings.clear();
186     names.clear();
187     for (dim_t i=0; i<m_mpiInfo->size; i++) {
188     stringstream path;
189     path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
190     tempstrings.push_back(path.str());
191     names.push_back((char*)tempstrings.back().c_str());
192     }
193     types.assign(m_mpiInfo->size, DB_QUADVAR);
194     DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
195     &types[0], NULL);
196     }
197    
198     if (m_mpiInfo->size > 1) {
199     #ifdef ESYS_MPI
200     PMPIO_HandOffBaton(baton, dbfile);
201     PMPIO_Finish(baton);
202     #endif
203     } else {
204     DBClose(dbfile);
205     }
206    
207     #else // USE_SILO
208     throw RipleyException("dump(): no Silo support");
209     #endif
210     }
211    
212     const int* Brick::borrowSampleReferenceIDs(int fsType) const
213     {
214 caltinay 3697 switch (fsType) {
215     case Nodes:
216     return &m_nodeId[0];
217     case Elements:
218     return &m_elementId[0];
219     case FaceElements:
220     return &m_faceId[0];
221     default:
222     break;
223     }
224 caltinay 3691
225 caltinay 3697 stringstream msg;
226     msg << "borrowSampleReferenceIDs() not implemented for function space type "
227     << fsType;
228     throw RipleyException(msg.str());
229 caltinay 3691 }
230    
231     bool Brick::ownSample(int fsCode, index_t id) const
232     {
233     #ifdef ESYS_MPI
234     if (fsCode == Nodes) {
235     const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
236     const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
237     return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
238     } else
239     throw RipleyException("ownSample() only implemented for Nodes");
240     #else
241     return true;
242     #endif
243     }
244    
245     void Brick::interpolateOnDomain(escript::Data& target,
246     const escript::Data& in) const
247     {
248     const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));
249     const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));
250     if (inDomain != *this)
251     throw RipleyException("Illegal domain of interpolant");
252     if (targetDomain != *this)
253     throw RipleyException("Illegal domain of interpolation target");
254    
255     throw RipleyException("interpolateOnDomain() not implemented");
256     }
257    
258     Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
259     bool reducedColOrder) const
260     {
261     if (reducedRowOrder || reducedColOrder)
262     throw RipleyException("getPattern() not implemented for reduced order");
263    
264     throw RipleyException("getPattern() not implemented");
265     }
266    
267     void Brick::Print_Mesh_Info(const bool full) const
268     {
269     RipleyDomain::Print_Mesh_Info(full);
270     if (full) {
271     cout << " Id Coordinates" << endl;
272     cout.precision(15);
273     cout.setf(ios::scientific, ios::floatfield);
274     for (index_t i=0; i < getNumNodes(); i++) {
275     cout << " " << setw(5) << m_nodeId[i]
276     << " " << (m_l0*(i%m_N0+m_offset0))/m_gNE0
277     << " " << (m_l1*(i%(m_N0*m_N1)/m_N0+m_offset1))/m_gNE1
278     << " " << (m_l2*(i/(m_N0*m_N1)+m_offset2))/m_gNE2 << endl;
279     }
280     }
281     }
282    
283     //protected
284     dim_t Brick::getNumFaceElements() const
285     {
286     dim_t n=0;
287     //left
288     if (m_offset0==0)
289     n+=m_NE1*m_NE2;
290     //right
291     if (m_mpiInfo->rank%m_NX==m_NX-1)
292     n+=m_NE1*m_NE2;
293     //bottom
294     if (m_offset1==0)
295     n+=m_NE0*m_NE2;
296     //top
297     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
298     n+=m_NE0*m_NE2;
299     //front
300     if (m_offset2==0)
301     n+=m_NE0*m_NE1;
302     //back
303     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
304     n+=m_NE0*m_NE1;
305    
306     return n;
307     }
308    
309     //protected
310     void Brick::assembleCoordinates(escript::Data& arg) const
311     {
312     escriptDataC x = arg.getDataC();
313     int numDim = m_numDim;
314     if (!isDataPointShapeEqual(&x, 1, &numDim))
315     throw RipleyException("setToX: Invalid Data object shape");
316     if (!numSamplesEqual(&x, 1, getNumNodes()))
317     throw RipleyException("setToX: Illegal number of samples in Data object");
318    
319     arg.requireWrite();
320     #pragma omp parallel for
321     for (dim_t i2 = 0; i2 < m_N2; i2++) {
322     for (dim_t i1 = 0; i1 < m_N1; i1++) {
323     for (dim_t i0 = 0; i0 < m_N0; i0++) {
324     double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
325     point[0] = (m_l0*(i0+m_offset0))/m_gNE0;
326     point[1] = (m_l1*(i1+m_offset1))/m_gNE1;
327     point[2] = (m_l2*(i2+m_offset2))/m_gNE2;
328     }
329     }
330     }
331     }
332    
333     //private
334     void Brick::populateSampleIds()
335     {
336 caltinay 3697 // identifiers are ordered from left to right, bottom to top, front to back
337     // on each rank, except for the shared nodes which are owned by the rank
338     // below / to the left / to the front of the current rank
339 caltinay 3691 const index_t firstId = getNumNodes()*m_mpiInfo->rank;
340     const index_t diff0 = m_N0*(m_N1*m_N2-1)+1;
341     const index_t diff1 = m_N0*m_N1*(m_N2*m_NX-1)+m_N0;
342     const index_t diff2 = m_N0*m_N1*m_N2*m_NX*m_NY-m_N0*m_N1*(m_N2-1);
343     m_nodeId.resize(getNumNodes());
344     #pragma omp parallel for
345     for (dim_t k=0; k<getNumNodes(); k++) {
346     index_t id = firstId+k;
347     if (m_offset0 > 0 && k%m_N0==0)
348     id -= diff0;
349     if (m_offset1 > 0 && k%(m_N0*m_N1)<m_N0)
350     id -= diff1;
351     if (m_offset2 > 0 && k/(m_N0*m_N1)==0)
352     id -= diff2;
353     m_nodeId[k]=id;
354     }
355 caltinay 3697
356     // elements
357     m_elementId.resize(getNumElements());
358     #pragma omp parallel for
359     for (dim_t k=0; k<getNumElements(); k++) {
360     m_elementId[k]=k;
361     }
362    
363     // face elements
364     m_faceId.resize(getNumFaceElements());
365     #pragma omp parallel for
366     for (dim_t k=0; k<getNumFaceElements(); k++) {
367     m_faceId[k]=k;
368     }
369 caltinay 3691 }
370    
371     } // end of namespace ripley
372    

  ViewVC Help
Powered by ViewVC 1.1.26