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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3691 - (show annotations)
Wed Nov 23 23:07:37 2011 UTC (7 years, 10 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
File size: 10904 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/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 #pragma omp parallel
142 {
143 #pragma omp for
144 for (dim_t i0 = 0; i0 < m_N0; i0++) {
145 coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;
146 }
147 #pragma omp for
148 for (dim_t i1 = 0; i1 < m_N1; i1++) {
149 coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;
150 }
151 }
152 int dims[] = { m_N0, m_N1 };
153 DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,
154 DB_COLLINEAR, NULL);
155
156 if (m_mpiInfo->rank == 0) {
157 vector<string> tempstrings;
158 vector<char*> meshnames;
159 for (dim_t i=0; i<m_mpiInfo->size; i++) {
160 stringstream path;
161 path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
162 tempstrings.push_back(path.str());
163 meshnames.push_back((char*)tempstrings.back().c_str());
164 }
165 vector<int> meshtypes(m_mpiInfo->size, DB_QUAD_RECT);
166 DBSetDir(dbfile, "/");
167 DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &meshnames[0],
168 &meshtypes[0], NULL);
169 }
170
171 if (m_mpiInfo->size > 1) {
172 #ifdef ESYS_MPI
173 PMPIO_HandOffBaton(baton, dbfile);
174 PMPIO_Finish(baton);
175 #endif
176 } else {
177 DBClose(dbfile);
178 }
179
180 #else // USE_SILO
181 throw RipleyException("dump(): no Silo support");
182 #endif
183 }
184
185 const int* Rectangle::borrowSampleReferenceIDs(int functionSpaceType) const
186 {
187 switch (functionSpaceType) {
188 case Nodes:
189 return &m_nodeId[0];
190 case Elements:
191 return &m_elementId[0];
192 case FaceElements:
193 return &m_faceId[0];
194 default:
195 break;
196 }
197
198 stringstream msg;
199 msg << "borrowSampleReferenceIDs() not implemented for function space type "
200 << functionSpaceType;
201 throw RipleyException(msg.str());
202 }
203
204 bool Rectangle::ownSample(int fsCode, index_t id) const
205 {
206 #ifdef ESYS_MPI
207 if (fsCode == Nodes) {
208 const index_t myFirst=getNumNodes()*m_mpiInfo->rank;
209 const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;
210 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
211 } else
212 throw RipleyException("ownSample() only implemented for Nodes");
213 #else
214 return true;
215 #endif
216 }
217
218 void Rectangle::interpolateOnDomain(escript::Data& target,
219 const escript::Data& in) const
220 {
221 const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));
222 const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));
223 if (inDomain != *this)
224 throw RipleyException("Illegal domain of interpolant");
225 if (targetDomain != *this)
226 throw RipleyException("Illegal domain of interpolation target");
227
228 throw RipleyException("interpolateOnDomain() not implemented");
229 }
230
231 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
232 bool reducedColOrder) const
233 {
234 if (reducedRowOrder || reducedColOrder)
235 throw RipleyException("getPattern() not implemented for reduced order");
236
237 /*
238 // create distribution
239 IndexVector dist;
240 for (index_t i=0; i<m_mpiInfo->size+1; i++)
241 dist.push_back(i*getNumNodes());
242 Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
243 &dist[0], 1, 0);
244
245 // connectors
246 dim_t numNeighbours = ...;
247 RankVector neighbour(numNeighbours);
248 IndexVector offsetInShared(numNeighbours+1);
249 IndexVector shared(offsetInShared[numNeighbours]);
250
251 Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
252 getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
253 1, 0, m_mpiInfo);
254 Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
255 getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,
256 1, 0, m_mpiInfo);
257 Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
258
259 // create patterns
260 Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
261 M, N, ptr, index);
262 Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
263 M, N, ...);
264 Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
265 ...);
266
267 Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
268 MATRIX_FORMAT_DEFAULT, distribution, distribution,
269 mainPattern, colCouplePattern, rowCouplePattern,
270 connector, connector);
271 Paso_Pattern_free(mainPattern);
272 Paso_Pattern_free(colCouplePattern);
273 Paso_Pattern_free(rowCouplePattern);
274 Paso_Distribution_free(distribution);
275 Paso_SharedComponents_free(snd_shcomp);
276 Paso_SharedComponents_free(rcv_shcomp);
277 */
278 throw RipleyException("getPattern() not implemented");
279 }
280
281 void Rectangle::Print_Mesh_Info(const bool full) const
282 {
283 RipleyDomain::Print_Mesh_Info(full);
284 if (full) {
285 cout << " Id Coordinates" << endl;
286 cout.precision(15);
287 cout.setf(ios::scientific, ios::floatfield);
288 for (index_t i=0; i < getNumNodes(); i++) {
289 cout << " " << setw(5) << m_nodeId[i]
290 << " " << (m_l0*(i%m_N0+m_offset0))/m_gNE0
291 << " " << (m_l1*(i/m_N0+m_offset1))/m_gNE1 << endl;
292 }
293 }
294 }
295
296 //protected
297 dim_t Rectangle::getNumFaceElements() const
298 {
299 dim_t n=0;
300 //left
301 if (m_offset0==0)
302 n+=m_NE1;
303 //right
304 if (m_mpiInfo->rank%m_NX==m_NX-1)
305 n+=m_NE1;
306 //bottom
307 if (m_offset1==0)
308 n+=m_NE0;
309 //top
310 if (m_mpiInfo->rank/m_NX==m_NY-1)
311 n+=m_NE0;
312
313 return n;
314 }
315
316 //protected
317 void Rectangle::assembleCoordinates(escript::Data& arg) const
318 {
319 escriptDataC x = arg.getDataC();
320 int numDim = m_numDim;
321 if (!isDataPointShapeEqual(&x, 1, &numDim))
322 throw RipleyException("setToX: Invalid Data object shape");
323 if (!numSamplesEqual(&x, 1, getNumNodes()))
324 throw RipleyException("setToX: Illegal number of samples in Data object");
325
326 arg.requireWrite();
327 #pragma omp parallel for
328 for (dim_t i1 = 0; i1 < m_N1; i1++) {
329 for (dim_t i0 = 0; i0 < m_N0; i0++) {
330 double* point = arg.getSampleDataRW(i0+m_N0*i1);
331 point[0] = (m_l0*(i0+m_offset0))/m_gNE0;
332 point[1] = (m_l1*(i1+m_offset1))/m_gNE1;
333 }
334 }
335 }
336
337 //private
338 void Rectangle::populateSampleIds()
339 {
340 const index_t firstId = getNumNodes()*m_mpiInfo->rank;
341 const index_t diff0 = m_N0*(m_N1-1)+1;
342 const index_t diff1 = m_N0*((m_NX-1)*m_N1+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==0)
350 id -= diff1;
351 m_nodeId[k]=id;
352 }
353 }
354
355 } // end of namespace ripley
356

  ViewVC Help
Powered by ViewVC 1.1.26