# Contents of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

Revision 3691 - (show annotations)
Wed Nov 23 23:07:37 2011 UTC (8 years, 3 months ago) by caltinay
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 7 * 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 11 * 12 *******************************************************/ 13 14 #include 15 extern "C" { 16 #include "paso/SystemMatrixPattern.h" 17 } 18 19 #if USE_SILO 20 #include 21 #ifdef ESYS_MPI 22 #include 23 #endif 24 #endif 25 26 #include 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(&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 x(new double[m_N0]); 139 boost::scoped_ptr 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 tempstrings; 158 vector meshnames; 159 for (dim_t i=0; isize; 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 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(*(in.getFunctionSpace().getDomain())); 222 const Rectangle& targetDomain=dynamic_cast(*(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; isize+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 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