/[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 3691 - (hide annotations)
Wed Nov 23 23:07:37 2011 UTC (7 years, 4 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 10411 byte(s)
Next iteration

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     DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 3, NULL, 0,
170     DB_INT, DB_NODECENT, NULL);
171    
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     if (fsType == Nodes)
215     return &m_nodeId[0];
216    
217     throw RipleyException("borrowSampleReferenceIDs() only implemented for Nodes");
218     }
219    
220     bool Brick::ownSample(int fsCode, index_t id) const
221     {
222     #ifdef ESYS_MPI
223     if (fsCode == Nodes) {
224     const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
225     const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
226     return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
227     } else
228     throw RipleyException("ownSample() only implemented for Nodes");
229     #else
230     return true;
231     #endif
232     }
233    
234     void Brick::interpolateOnDomain(escript::Data& target,
235     const escript::Data& in) const
236     {
237     const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));
238     const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));
239     if (inDomain != *this)
240     throw RipleyException("Illegal domain of interpolant");
241     if (targetDomain != *this)
242     throw RipleyException("Illegal domain of interpolation target");
243    
244     throw RipleyException("interpolateOnDomain() not implemented");
245     }
246    
247     Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
248     bool reducedColOrder) const
249     {
250     if (reducedRowOrder || reducedColOrder)
251     throw RipleyException("getPattern() not implemented for reduced order");
252    
253     throw RipleyException("getPattern() not implemented");
254     }
255    
256     void Brick::Print_Mesh_Info(const bool full) const
257     {
258     RipleyDomain::Print_Mesh_Info(full);
259     if (full) {
260     cout << " Id Coordinates" << endl;
261     cout.precision(15);
262     cout.setf(ios::scientific, ios::floatfield);
263     for (index_t i=0; i < getNumNodes(); i++) {
264     cout << " " << setw(5) << m_nodeId[i]
265     << " " << (m_l0*(i%m_N0+m_offset0))/m_gNE0
266     << " " << (m_l1*(i%(m_N0*m_N1)/m_N0+m_offset1))/m_gNE1
267     << " " << (m_l2*(i/(m_N0*m_N1)+m_offset2))/m_gNE2 << endl;
268     }
269     }
270     }
271    
272     //protected
273     dim_t Brick::getNumFaceElements() const
274     {
275     dim_t n=0;
276     //left
277     if (m_offset0==0)
278     n+=m_NE1*m_NE2;
279     //right
280     if (m_mpiInfo->rank%m_NX==m_NX-1)
281     n+=m_NE1*m_NE2;
282     //bottom
283     if (m_offset1==0)
284     n+=m_NE0*m_NE2;
285     //top
286     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
287     n+=m_NE0*m_NE2;
288     //front
289     if (m_offset2==0)
290     n+=m_NE0*m_NE1;
291     //back
292     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
293     n+=m_NE0*m_NE1;
294    
295     return n;
296     }
297    
298     //protected
299     void Brick::assembleCoordinates(escript::Data& arg) const
300     {
301     escriptDataC x = arg.getDataC();
302     int numDim = m_numDim;
303     if (!isDataPointShapeEqual(&x, 1, &numDim))
304     throw RipleyException("setToX: Invalid Data object shape");
305     if (!numSamplesEqual(&x, 1, getNumNodes()))
306     throw RipleyException("setToX: Illegal number of samples in Data object");
307    
308     arg.requireWrite();
309     #pragma omp parallel for
310     for (dim_t i2 = 0; i2 < m_N2; i2++) {
311     for (dim_t i1 = 0; i1 < m_N1; i1++) {
312     for (dim_t i0 = 0; i0 < m_N0; i0++) {
313     double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
314     point[0] = (m_l0*(i0+m_offset0))/m_gNE0;
315     point[1] = (m_l1*(i1+m_offset1))/m_gNE1;
316     point[2] = (m_l2*(i2+m_offset2))/m_gNE2;
317     }
318     }
319     }
320     }
321    
322     //private
323     void Brick::populateSampleIds()
324     {
325     const index_t firstId = getNumNodes()*m_mpiInfo->rank;
326     const index_t diff0 = m_N0*(m_N1*m_N2-1)+1;
327     const index_t diff1 = m_N0*m_N1*(m_N2*m_NX-1)+m_N0;
328     const index_t diff2 = m_N0*m_N1*m_N2*m_NX*m_NY-m_N0*m_N1*(m_N2-1);
329     m_nodeId.resize(getNumNodes());
330     #pragma omp parallel for
331     for (dim_t k=0; k<getNumNodes(); k++) {
332     index_t id = firstId+k;
333     if (m_offset0 > 0 && k%m_N0==0)
334     id -= diff0;
335     if (m_offset1 > 0 && k%(m_N0*m_N1)<m_N0)
336     id -= diff1;
337     if (m_offset2 > 0 && k/(m_N0*m_N1)==0)
338     id -= diff2;
339     m_nodeId[k]=id;
340     }
341     }
342    
343     } // end of namespace ripley
344    

  ViewVC Help
Powered by ViewVC 1.1.26