- Timestamp:
- 06/06/17 17:58:16 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/dev_olga/src/transformation/domain_algorithm_expand.cpp
r978 r1158 45 45 CDomain* domainSource, 46 46 CExpandDomain* expandDomain) 47 : CDomainAlgorithmTransformation(domainDestination, domainSource) 47 : CDomainAlgorithmTransformation(domainDestination, domainSource), 48 isXPeriodic_(false), isYPeriodic_(false) 48 49 { 49 50 if (domainDestination == domainSource) … … 55 56 } 56 57 57 // if (!domainDestination->hasRefTo(domainSource))58 // {59 // ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",60 // << "Domain domain destination must refer to domain source (directly or indirectly) by domain_ref " << std::endl61 // << "Domain source " <<domainSource->getId() << std::endl62 // << "Domain destination " <<domainDestination->getId() << std::endl);63 // }64 65 58 this->type_ = (ELEMENT_MODIFICATION_WITH_DATA); 59 // Make sure domain source have all valid attributes 60 // domainSource->checkAllAttributes(); 66 61 expandDomain->checkValid(domainDestination); 67 62 if (!expandDomain->i_periodic.isEmpty()) isXPeriodic_ = expandDomain->i_periodic; 63 if (!expandDomain->j_periodic.isEmpty()) isYPeriodic_ = expandDomain->j_periodic; 64 68 65 switch (expandDomain->type) 69 66 { … … 92 89 CContextClient* client=context->client; 93 90 91 int type = 1; // For edge 92 CMesh mesh; 94 93 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 95 94 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 96 95 CArray<int,2> neighborsSrc; 97 98 int type = 1; // For edge 99 CMesh mesh; 100 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 101 updateDomainAttributes(domainDestination, domainSource, neighborsSrc); 96 switch (domainSource->type) { 97 case CDomain::type_attr::unstructured: 98 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 99 updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc); 100 break; 101 default: 102 updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc); 103 break; 104 } 102 105 } 103 106 … … 113 116 CContextClient* client=context->client; 114 117 118 int type = 1; // For edge 119 CMesh mesh; 115 120 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 116 121 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 117 122 CArray<int,2> neighborsSrc; 118 119 int type = 0; // For node 120 CMesh mesh; 121 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 122 updateDomainAttributes(domainDestination, domainSource, neighborsSrc); 123 switch (domainSource->type) { 124 case CDomain::type_attr::unstructured: 125 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 126 updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc); 127 break; 128 default: 129 updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc); 130 break; 131 } 132 } 133 134 /*! 135 * Extend rectilinear or curvilinear domain destination and update its attributes 136 * Suppose that domain destination and domain source have the same values for all attributes (by inheritance) 137 * \param [in/out] domainDestination domain destination 138 * \param [in] domainSource domain source 139 * \param [in] neighborsDomainSrc neighbor of domain source. For now, we don't need it for rectilinear 140 */ 141 void CDomainAlgorithmExpand::updateRectilinearDomainAttributes(CDomain* domainDestination, 142 CDomain* domainSource, 143 CArray<int,2>& neighborsDomainSrc) 144 { 145 int index, globalIndex, idx; 146 int iindexDst, jindexDst, globIndexDst; 147 int iindexSrc, jindexSrc, globIndexSrc; 148 CContext* context = CContext::getCurrent(); 149 CContextClient* client=context->client; 150 151 // First of all, "copy" all attributes of domain source to domain destination 152 StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue() 153 : ""; 154 if (domainDstRef != domainSource->getId()) 155 { 156 domainDestination->domain_ref.setValue(domainSource->getId()); 157 domainDestination->solveRefInheritance(true); 158 } 159 160 if (domainDstRef.empty()) domainDestination->domain_ref.reset(); 161 else domainDestination->domain_ref.setValue(domainDstRef); 162 163 // Here are attributes of source need tranfering 164 int niGloSrc = domainSource->ni_glo; 165 int njGloSrc = domainSource->nj_glo; 166 int niSrc = domainSource->ni, ibegin = domainSource->ibegin; 167 int njSrc = domainSource->nj, jbegin = domainSource->jbegin; 168 int dataDimSrc = domainSource->data_dim; 169 CArray<bool,1>& mask_1d_src = domainSource->mask_1d; 170 CArray<int,1>& i_index_src = domainSource->i_index; 171 CArray<int,1>& j_index_src = domainSource->j_index; 172 CArray<int,1>& data_i_index_src = domainSource->data_i_index; 173 CArray<int,1>& data_j_index_src = domainSource->data_j_index; 174 int data_i_begin_src = domainSource->data_ibegin; 175 int data_j_begin_src = domainSource->data_jbegin; 176 CArray<double,1>& lon_src = domainSource->lonvalue; 177 CArray<double,1>& lat_src = domainSource->latvalue; 178 179 // We need to generate boundary for longitude and latitude 180 if (domainSource->bounds_lon_1d.isEmpty() || domainSource->bounds_lat_1d.isEmpty()) 181 { 182 CArray<double,1> lon = lon_src(Range(0,niSrc-1)); 183 CArray<double,1> lat = lat_src(Range(0,lat_src.numElements()-niSrc,niSrc)); 184 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 185 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 186 domainSource->fillInRectilinearBoundLonLat(lon_src, lat_src, bounds_lon_src, bounds_lat_src); 187 } 188 189 190 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 191 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 192 193 int nVertex = bounds_lon_src.shape()[0]; 194 int oldNbLocal = i_index_src.numElements(); 195 // Calculate ni, nj by using i_index and j_index 196 int niSrcByIndex = max(i_index_src) - min(i_index_src) + 1; 197 int njSrcByIndex = max(j_index_src) - min(j_index_src) + 1; 198 int dataIindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : niSrcByIndex; 199 int dataJindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : njSrcByIndex; 200 201 // Uncompress data_i_index, data_j_index 202 CArray<int,1> data_i_index_src_full(oldNbLocal); 203 CArray<int,1> data_j_index_src_full(oldNbLocal); 204 int nbUnMaskedPointOnLocalDomain = 0; 205 data_i_index_src_full = -1; // Suppose all values are masked 206 data_j_index_src_full = -1; // Suppose all values are masked 207 for (idx = 0; idx < data_i_index_src.numElements(); ++idx) 208 { 209 int dataIidx = data_i_index_src(idx) + data_i_begin_src; 210 int dataJidx = data_j_index_src(idx) + data_j_begin_src; 211 if ((0 <= dataIidx) && (dataIidx < dataIindexBoundSrc) && 212 (0 <= dataJidx) && (dataJidx < dataJindexBoundSrc)) 213 { 214 data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIidx; 215 data_j_index_src_full(nbUnMaskedPointOnLocalDomain) = dataJidx; 216 ++nbUnMaskedPointOnLocalDomain; 217 } 218 } 219 220 // Expand domain destination, not only local but also global 221 int niGloDst = niGloSrc + 2; 222 int njGloDst = njGloSrc + 2; 223 int niDst = niSrc + 2; 224 int njDst = njSrc + 2; 225 domainDestination->ni_glo.setValue(niGloDst); 226 domainDestination->nj_glo.setValue(njGloDst); 227 domainDestination->ni.setValue(niDst); 228 domainDestination->nj.setValue(njDst); 229 domainDestination->global_zoom_ni.setValue(domainSource->global_zoom_ni+2); 230 domainDestination->global_zoom_nj.setValue(domainSource->global_zoom_nj+2); 231 232 CArray<bool,1>& mask_1d_dst = domainDestination->mask_1d; 233 CArray<int,1>& i_index_dst = domainDestination->i_index; 234 CArray<int,1>& j_index_dst = domainDestination->j_index; 235 CArray<int,1>& data_i_index_dst = domainDestination->data_i_index; 236 CArray<int,1>& data_j_index_dst = domainDestination->data_j_index; 237 238 // Make sure that we use only lonvalue_client, latvalue_client 239 if (!domainDestination->lonvalue_1d.isEmpty()) domainDestination->lonvalue_1d.reset(); 240 if (!domainDestination->latvalue_1d.isEmpty()) domainDestination->latvalue_1d.reset(); 241 if (!domainDestination->lonvalue_2d.isEmpty()) domainDestination->lonvalue_2d.reset(); 242 if (!domainDestination->latvalue_2d.isEmpty()) domainDestination->latvalue_2d.reset(); 243 244 // Recalculate i_index, j_index of extended domain 245 // Should be enough for common case, but if we have arbitrary distribution? 246 int newNbLocalDst = niDst * njDst; 247 248 mask_1d_dst.resize(newNbLocalDst); 249 i_index_dst.resize(newNbLocalDst); 250 j_index_dst.resize(newNbLocalDst); 251 CArray<int,1> data_i_index_dst_full(newNbLocalDst); 252 CArray<int,1> data_j_index_dst_full(newNbLocalDst); 253 254 domainDestination->lonvalue.resizeAndPreserve(newNbLocalDst); 255 domainDestination->latvalue.resizeAndPreserve(newNbLocalDst); 256 domainDestination->bounds_lon_1d.resizeAndPreserve(nVertex, newNbLocalDst); 257 domainDestination->bounds_lat_1d.resizeAndPreserve(nVertex, newNbLocalDst); 258 CArray<double,1>& lon_dst = domainDestination->lonvalue; 259 CArray<double,1>& lat_dst = domainDestination->latvalue; 260 CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d; 261 CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d; 262 263 // Update i_index, j_index 264 for (int j = 0; j < njDst; ++j) 265 for (int i = 0; i < niDst; ++i) 266 { 267 idx = j * niDst + i; 268 i_index_dst(idx) = i + ibegin; 269 j_index_dst(idx) = j + jbegin; 270 } 271 272 273 // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...) 274 // Global index mapping between destination and source 275 this->transformationMapping_.resize(1); 276 this->transformationWeight_.resize(1); 277 TransformationIndexMap& transMap = this->transformationMapping_[0]; 278 TransformationWeightMap& transWeight = this->transformationWeight_[0]; 279 280 transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor())); 281 transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor())); 282 283 // Index mapping for local domain 284 // Mapping global index of expanded domain into original one 285 // (Representing global index of expanded domain in form of global index of original one) 286 CArray<size_t,1> globalIndexSrcOnDstDomain(newNbLocalDst); 287 for (idx = 0; idx < newNbLocalDst; ++idx) 288 { 289 iindexDst = i_index_dst(idx); 290 jindexDst = j_index_dst(idx); 291 globIndexDst = jindexDst * niGloDst + iindexDst; 292 globIndexSrc = (((jindexDst-1)+njGloSrc) % njGloSrc) * niGloSrc + (((iindexDst-1)+niGloSrc) % niGloSrc) ; 293 globalIndexSrcOnDstDomain(idx) = globIndexSrc; 294 295 transMap[globIndexDst].push_back(globIndexSrc); 296 transWeight[globIndexDst].push_back(1.0); 297 } 298 299 // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...) 300 CClientClientDHTDouble::Index2VectorInfoTypeMap localData; 301 localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor())); 302 303 // Information exchanged among domains (attention to their order), number in parentheses presents size of data 304 // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1) 305 int dataPackageSize = 1 + 1 + // lon + lat 306 nVertex + nVertex + //bounds_lon + bounds_lat 307 1 + // mask_1d_dst; 308 1 + 1; // data_i_index + data_j_index 309 // Initialize database 310 for (int idx = 0; idx < oldNbLocal; ++idx) 311 { 312 index = i_index_src(idx) + j_index_src(idx) * niGloSrc; 313 localData[index].resize(dataPackageSize); 314 std::vector<double>& data = localData[index]; 315 316 //Pack data 317 int dataIdx = 0; 318 data[dataIdx] = lon_src(idx);++dataIdx; 319 data[dataIdx] = lat_src(idx);++dataIdx; 320 for (int i = 0; i < nVertex; ++i) 321 { 322 data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx; 323 } 324 for (int i = 0; i < nVertex; ++i) 325 { 326 data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx; 327 } 328 data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx; 329 data[dataIdx] = data_i_index_src_full(idx);++dataIdx; 330 data[dataIdx] = data_j_index_src_full(idx); 331 } 332 333 CClientClientDHTDouble dhtData(localData,client->intraComm); 334 dhtData.computeIndexInfoMapping(globalIndexSrcOnDstDomain); 335 CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap(); 336 CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it; 337 338 // Ok get all data for destination 339 // If domain is not periodic, then we mask all extended part. 340 int nbUnMaskedPointOnExtendedPart = 0, remainder = 0, dataIIndex, dataJIndex; 341 size_t nIdx; 342 double maskValue = 1.0; 343 for (index = 0; index < newNbLocalDst; ++index) 344 { 345 nIdx = globalIndexSrcOnDstDomain(index); 346 it = neighborData.find(nIdx); 347 if (ite != it) 348 { 349 std::vector<double>& data = it->second; 350 // Unpack data 351 int dataIdx = 0; 352 lon_dst(index) = data[dataIdx]; ++dataIdx; 353 lat_dst(index) = data[dataIdx]; ++dataIdx; 354 for (int i = 0; i < nVertex; ++i) 355 { 356 bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx; 357 } 358 for (int i = 0; i < nVertex; ++i) 359 { 360 bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx; 361 } 362 363 // Check whether we have x periodic. If we don't, we should mask all point at 0 and niGloDst-1 364 maskValue = data[dataIdx]; 365 if (!isXPeriodic_) 366 { 367 remainder = i_index_dst(index) % (niGloDst-1); 368 if (0 == remainder) 369 { 370 maskValue = -1.0; 371 } 372 } 373 374 if (!isYPeriodic_) 375 { 376 remainder = j_index_dst(index) % (njGloDst-1); 377 if (0 == remainder) 378 { 379 maskValue = -1.0; 380 } 381 } 382 383 mask_1d_dst(index) = (1.0 == maskValue) ? true : false; ++dataIdx; 384 385 dataIIndex = (int) data[dataIdx]; 386 if (!isXPeriodic_) 387 { 388 remainder = i_index_dst(index) % (niGloDst-1); 389 if (0 == remainder) 390 { 391 dataIIndex = -1; 392 } 393 } 394 data_i_index_dst_full(index) = dataIIndex; ++dataIdx; 395 396 dataJIndex = (int) data[dataIdx]; 397 if (!isYPeriodic_) 398 { 399 remainder = j_index_dst(index) % (njGloDst-1); 400 if (0 == remainder) 401 { 402 dataJIndex = -1; 403 } 404 } 405 data_j_index_dst_full(index) = dataJIndex; 406 407 if ((0 <= data_i_index_dst_full(index)) && (0 <= data_j_index_dst_full(index))) 408 { 409 ++nbUnMaskedPointOnExtendedPart; 410 } 411 } 412 } 413 414 415 // Finally, update data_i_index, data_j_index 416 int dataDstDim = domainDestination->data_dim; 417 data_i_index_dst.resize(nbUnMaskedPointOnExtendedPart); 418 data_j_index_dst.resize(nbUnMaskedPointOnExtendedPart); 419 int count = 0; 420 for (idx = 0; idx < newNbLocalDst; ++idx) 421 { 422 dataIIndex = data_i_index_dst_full(idx); 423 dataJIndex = data_j_index_dst_full(idx); 424 if ((0 <= dataIIndex) && (0 <= dataJIndex)) 425 { 426 data_i_index_dst(count) = (1 == dataDstDim) ? idx : i_index_dst(idx) - i_index_dst(0); 427 data_j_index_dst(count) = (1 == dataDstDim) ? 0 : j_index_dst(idx) - j_index_dst(0); 428 ++count; 429 } 430 } 431 432 // Update data_ni, data_nj 433 434 domainDestination->data_ni.setValue((1==dataDstDim) ? niDst * njDst : niDst); 435 domainDestination->data_nj.setValue((1==dataDstDim) ? niDst * njDst : njDst); 436 domainDestination->data_ibegin.setValue(0); 437 domainDestination->data_jbegin.setValue(0); 438 439 // Update longitude and latitude 440 if (niSrc == domainSource->lonvalue_1d.numElements() && njSrc == domainSource->latvalue_1d.numElements()) // Ok, we have rectilinear here 441 { 442 domainDestination->lonvalue_1d.resize(niDst); 443 domainDestination->lonvalue_1d = lon_dst(Range(0,niDst-1)); 444 domainDestination->latvalue_1d.resize(njDst); 445 domainDestination->latvalue_1d = lat_dst(Range(0,lat_dst.numElements()-niDst,niDst)); 446 } 447 else // It should be curvilinear 448 { 449 domainDestination->lonvalue_1d.resize(lon_dst.numElements()); 450 domainDestination->lonvalue_1d = lon_dst; 451 domainDestination->latvalue_1d.resize(lat_dst.numElements()); 452 domainDestination->latvalue_1d = (lat_dst); 453 } 454 123 455 } 124 456 … … 130 462 * \param [in] neighborsDomainSrc domain extended part 131 463 */ 132 void CDomainAlgorithmExpand::updateDomainAttributes(CDomain* domainDestination, 133 CDomain* domainSource, 134 CArray<int,2>& neighborsDomainSrc) 135 { 464 void CDomainAlgorithmExpand::updateUnstructuredDomainAttributes(CDomain* domainDestination, 465 CDomain* domainSource, 466 CArray<int,2>& neighborsDomainSrc) 467 { 468 136 469 CContext* context = CContext::getCurrent(); 137 470 CContextClient* client=context->client; … … 255 588 data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx; 256 589 } 257 data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1 ; ++dataIdx;590 data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx; 258 591 data[dataIdx] = data_i_index_src_full(idx); 259 592 } … … 308 641 { 309 642 dataIdx = data_i_index_dst_full(idx); 310 if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))643 if ((0 <= dataIdx)) 311 644 { 312 645 ++count; … … 322 655 { 323 656 dataIdx = data_i_index_dst_full(idx); 324 if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))657 if ((0 <= dataIdx)) 325 658 { 326 659 data_i_index_dst(count) = dataIdx; … … 331 664 // Update ni 332 665 domainDestination->ni.setValue(newNbLocalDst); 333 334 666 } 335 667
Note: See TracChangeset
for help on using the changeset viewer.