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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3691 - (hide annotations)
Wed Nov 23 23:07:37 2011 UTC (7 years, 11 months ago) by caltinay
File size: 10904 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/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