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

Contents of /trunk/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3691 - (show annotations)
Wed Nov 23 23:07:37 2011 UTC (7 years, 5 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 10411 byte(s)
Next iteration

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/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