Changeset 2011 for XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/generic_algorithm_transformation.cpp
- Timestamp:
- 01/12/21 23:05:02 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/generic_algorithm_transformation.cpp
r2007 r2011 23 23 24 24 CGenericAlgorithmTransformation::CGenericAlgorithmTransformation(bool isSource) 25 : isSource_(isSource), transformationMapping_(), transformationWeight_(), transformationPosition_(), 26 idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA), indexElementSrc_(), 27 computedProcSrcNonTransformedElement_(false), eliminateRedondantSrc_(true), isDistributedComputed_(false) 25 : isSource_(isSource) 28 26 { 29 27 } 30 28 31 void CGenericAlgorithmTransformation::updateData(CArray<double,1>& dataOut)32 {33 34 }35 36 void CGenericAlgorithmTransformation::apply(const std::vector<std::pair<int,double> >& localIndex,37 const double* dataInput,38 CArray<double,1>& dataOut,39 std::vector<bool>& flagInitial,40 bool ignoreMissingValue, bool firstPass )41 TRY42 {43 int nbLocalIndex = localIndex.size();44 double defaultValue = std::numeric_limits<double>::quiet_NaN();45 46 if (ignoreMissingValue)47 {48 if (firstPass) dataOut=defaultValue ;49 50 for (int idx = 0; idx < nbLocalIndex; ++idx)51 {52 if (! NumTraits<double>::isNan(*(dataInput + idx)))53 {54 if (flagInitial[localIndex[idx].first]) dataOut(localIndex[idx].first) = *(dataInput + idx) * localIndex[idx].second;55 else dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;56 flagInitial[localIndex[idx].first] = false; // Reset flag to indicate not all data source are nan57 }58 }59 60 }61 else62 {63 for (int idx = 0; idx < nbLocalIndex; ++idx)64 {65 dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;66 }67 }68 }69 CATCH70 71 void CGenericAlgorithmTransformation::computePositionElements(CGrid* dst, CGrid* src)72 TRY73 {74 int idxScalar = 0, idxAxis = 0, idxDomain = 0;75 CArray<int,1> axisDomainOrderDst = dst->axis_domain_order;76 for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)77 {78 int dimElement = axisDomainOrderDst(i);79 if (2 == dimElement)80 {81 elementPositionInGridDst2DomainPosition_[i] = idxDomain;82 ++idxDomain;83 }84 else if (1 == dimElement)85 {86 elementPositionInGridDst2AxisPosition_[i] = idxAxis;87 ++idxAxis;88 }89 else90 {91 elementPositionInGridDst2ScalarPosition_[i] = idxScalar;92 ++idxScalar;93 }94 }95 96 idxScalar = idxAxis = idxDomain = 0;97 CArray<int,1> axisDomainOrderSrc = src->axis_domain_order;98 for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i)99 {100 int dimElement = axisDomainOrderSrc(i);101 if (2 == dimElement)102 {103 elementPositionInGridSrc2DomainPosition_[i] = idxDomain;104 ++idxDomain;105 }106 else if (1 == dimElement)107 {108 elementPositionInGridSrc2AxisPosition_[i] = idxAxis;109 ++idxAxis;110 }111 else112 {113 elementPositionInGridSrc2ScalarPosition_[i] = idxScalar;114 ++idxScalar;115 }116 }117 }118 CATCH119 120 bool CGenericAlgorithmTransformation::isDistributedTransformation(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst)121 TRY122 {123 124 if (!isDistributedComputed_)125 {126 isDistributedComputed_=true ;127 if (!eliminateRedondantSrc_) isDistributed_=true ;128 else129 {130 CContext* context = CContext::getCurrent();131 132 computePositionElements(gridSrc, gridDst);133 std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();134 std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();135 std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();136 int distributed, distributed_glo ;137 138 CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;139 if (2 == axisDomainSrcOrder(elementPositionInGrid)) // It's domain140 {141 distributed=domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]->isDistributed() ;142 MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, context->intraComm_) ;143 144 }145 else if (1 == axisDomainSrcOrder(elementPositionInGrid))//it's an axis146 {147 distributed=axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]->isDistributed() ;148 MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, context->intraComm_) ;149 }150 else //it's a scalar151 {152 distributed_glo=false ;153 }154 isDistributed_=distributed_glo ;155 }156 }157 return isDistributed_ ;158 }159 CATCH160 161 /*!162 This function computes the global indexes of grid source, which the grid destination is in demand.163 \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order),164 then position of axis in grid is 1 and domain is positioned at 0.165 \param[in] gridSrc Grid source166 \param[in] gridDst Grid destination167 \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination168 */169 void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid,170 CGrid* gridSrc,171 CGrid* gridDst,172 SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)173 TRY174 {175 CContext* context = CContext::getCurrent();176 int nbClient = context->intraCommSize_;177 178 typedef std::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap;179 int idx;180 181 // compute position of elements on grids182 computePositionElements(gridDst, gridSrc);183 std::vector<CScalar*> scalarListDestP = gridDst->getScalars();184 std::vector<CAxis*> axisListDestP = gridDst->getAxis();185 std::vector<CDomain*> domainListDestP = gridDst->getDomains();186 CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;187 std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();188 std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();189 std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();190 CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;191 192 bool isTransPosEmpty = transformationPosition_.empty();193 CArray<size_t,1> transPos;194 if (!isTransPosEmpty) transPos.resize(transformationMapping_.size());195 std::set<size_t> allIndexSrc; // All index of source, which can be scattered among processes, need for doing transformation196 197 for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)198 {199 TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,200 iteTransMap = transformationMapping_[idxTrans].end();201 TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;202 203 // Build mapping between global source element index and global destination element index.204 itTransWeight = itbTransWeight;205 for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)206 {207 const std::vector<int>& srcIndex = itTransMap->second;208 for (idx = 0; idx < srcIndex.size(); ++idx)209 allIndexSrc.insert(srcIndex[idx]);210 }211 212 if (!isTransPosEmpty)213 {214 TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin();215 transPos(idxTrans) = itPosMap->second[0];216 }217 }218 219 size_t indexSrcSize = 0;220 CArray<size_t,1> indexSrc(allIndexSrc.size());221 for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it, ++indexSrcSize)222 indexSrc(indexSrcSize) = *it;223 224 // Flag to indicate whether we will recompute distribution of source global index on processes225 bool computeGlobalIndexOnProc = false;226 if (indexElementSrc_.size() != allIndexSrc.size())227 computeGlobalIndexOnProc = true;228 else229 {230 for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it)231 if (0 == indexElementSrc_.count(*it))232 {233 computeGlobalIndexOnProc = true;234 break;235 }236 }237 238 if (computeGlobalIndexOnProc)239 indexElementSrc_.swap(allIndexSrc);240 241 int sendValue = (computeGlobalIndexOnProc) ? 1 : 0;242 int recvValue = 0;243 MPI_Allreduce(&sendValue, &recvValue, 1, MPI_INT, MPI_SUM, context->intraComm_);244 computeGlobalIndexOnProc = (0 < recvValue);245 246 // CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc;247 248 if (computeGlobalIndexOnProc || !computedProcSrcNonTransformedElement_)249 {250 {251 CClientClientDHTInt::Index2VectorInfoTypeMap tmp ;252 globalIndexOfTransformedElementOnProc_.swap(tmp) ;253 }254 // Find out global index source of transformed element on corresponding process.255 if (globalElementIndexOnProc_.empty())256 globalElementIndexOnProc_.resize(axisDomainDstOrder.numElements());257 258 for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)259 {260 261 if (idx == elementPositionInGrid)262 computeExchangeGlobalIndex(indexSrc, axisDomainSrcOrder(idx), globalIndexOfTransformedElementOnProc_); //globalElementIndexOnProc[idx]);263 if (!computedProcSrcNonTransformedElement_)264 {265 if (2 == axisDomainDstOrder(idx)) // It's domain266 {267 if (idx != elementPositionInGrid)268 computeExchangeDomainIndex(domainListDestP[elementPositionInGridDst2DomainPosition_[idx]],269 domainListSrcP[elementPositionInGridSrc2DomainPosition_[idx]],270 transPos,271 globalElementIndexOnProc_[idx]);272 273 }274 else if (1 == axisDomainDstOrder(idx))//it's an axis275 {276 if (idx != elementPositionInGrid)277 computeExchangeAxisIndex(axisListDestP[elementPositionInGridDst2AxisPosition_[idx]],278 axisListSrcP[elementPositionInGridSrc2AxisPosition_[idx]],279 transPos,280 globalElementIndexOnProc_[idx]);281 }282 else //it's a scalar283 {284 if (idx != elementPositionInGrid)285 computeExchangeScalarIndex(scalarListDestP[elementPositionInGridDst2ScalarPosition_[idx]],286 scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[idx]],287 transPos,288 globalElementIndexOnProc_[idx]);289 290 }291 }292 }293 294 if (!isTransPosEmpty && !computedProcSrcNonTransformedElement_)295 {296 for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)297 {298 if (idx != elementPositionInGrid)299 {300 std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,301 ite = globalElementIndexOnProc_[idx].end();302 for (it = itb; it != ite; ++it) it->second.resize(1);303 }304 }305 }306 307 /*308 if (!computedProcSrcNonTransformedElement_)309 {310 for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)311 {312 if (idx != elementPositionInGrid)313 {314 std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,315 ite = globalElementIndexOnProc_[idx].end();316 for (it = itb; it != ite; ++it) procOfNonTransformedElements_.insert(it->first);317 if (procOfNonTransformedElements_.size() == nbClient)318 break;319 }320 }321 }322 323 // Processes contain the source index of transformed element324 std::set<int> procOfTransformedElement;325 CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc.begin(),326 itIdxe = globalIndexOfTransformedElementOnProc.end(), itIdx;327 for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx)328 {329 std::vector<int>& tmp = itIdx->second;330 for (int i = 0; i < tmp.size(); ++i)331 procOfTransformedElement.insert(tmp[i]);332 if (tmp.size() == nbClient)333 break;334 }335 336 std::set<int>& commonProc = (procOfTransformedElement.size() < procOfNonTransformedElements_.size()) ? procOfTransformedElement337 : (!procOfNonTransformedElements_.empty() ? procOfNonTransformedElements_ : procOfTransformedElement);338 339 std::vector<int> procContainSrcElementIdx(commonProc.size());340 int count = 0;341 for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it)342 procContainSrcElementIdx[count++] = *it;343 344 procContainSrcElementIdx_.swap(procContainSrcElementIdx);345 */346 347 if (procElementList_.empty()) procElementList_.resize(axisDomainDstOrder.numElements()) ;348 for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)349 {350 std::set<int>& procList=procElementList_[idx] ;351 std::set<int> commonTmp ;352 if (idx == elementPositionInGrid)353 {354 set<int> tmpSet ;355 procList.swap(tmpSet) ;356 CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc_.begin(),357 itIdxe = globalIndexOfTransformedElementOnProc_.end(), itIdx;358 for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx)359 {360 std::vector<int>& tmp = itIdx->second;361 for (int i = 0; i < tmp.size(); ++i) procList.insert(tmp[i]);362 if (tmp.size() == nbClient)363 break;364 }365 }366 else367 {368 if (!computedProcSrcNonTransformedElement_)369 {370 set<int> tmpSet ;371 procList.swap(tmpSet) ;372 std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,373 ite = globalElementIndexOnProc_[idx].end();374 for (it = itb; it != ite; ++it)375 {376 procList.insert(it->first);377 if (procList.size() == nbClient) break;378 }379 }380 }381 382 if (idx==0) commonProc_= procList ;383 else384 {385 for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it)386 if (procList.count(*it)==1) commonTmp.insert(*it) ;387 commonProc_.swap(commonTmp) ;388 }389 }390 std::vector<int> procContainSrcElementIdx(commonProc_.size());391 int count = 0;392 for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it) procContainSrcElementIdx[count++] = *it;393 procContainSrcElementIdx_.swap(procContainSrcElementIdx);394 395 // For the first time, surely we calculate proc containing non transformed elements396 if (!computedProcSrcNonTransformedElement_) computedProcSrcNonTransformedElement_ = true;397 }398 399 for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)400 {401 TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,402 iteTransMap = transformationMapping_[idxTrans].end();403 TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;404 SrcToDstMap src2DstMap;405 src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor()));406 407 // Build mapping between global source element index and global destination element index.408 std::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc_[elementPositionInGrid]);409 std::set<int> tmpCounter;410 itTransWeight = itbTransWeight;411 for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)412 {413 const std::vector<int>& srcIndex = itTransMap->second;414 const std::vector<double>& weight = itTransWeight->second;415 for (idx = 0; idx < srcIndex.size(); ++idx)416 {417 src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx]));418 if (0 == tmpCounter.count(srcIndex[idx]))419 {420 tmpCounter.insert(srcIndex[idx]);421 422 vector<int>& rankSrc = globalIndexOfTransformedElementOnProc_[srcIndex[idx]] ;423 for (int n=0;n<rankSrc.size();++n)424 {425 if (commonProc_.count(rankSrc[n])==1) globalElementIndexOnProc_[elementPositionInGrid][rankSrc[n]].push_back(srcIndex[idx]);426 }427 // for (int j = 0; j < procContainSrcElementIdx_.size(); ++j)428 // globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]);429 }430 }431 }432 433 if (!isTransPosEmpty)434 {435 for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)436 {437 if (idx != elementPositionInGrid)438 {439 std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,440 ite = globalElementIndexOnProc_[idx].end();441 for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans);442 }443 }444 }445 446 // Ok, now compute global index of grid source and ones of grid destination447 computeGlobalGridIndexMapping(elementPositionInGrid,448 procContainSrcElementIdx_, //srcRank,449 src2DstMap,450 gridSrc,451 gridDst,452 globalElementIndexOnProc_,453 globaIndexWeightFromSrcToDst);454 }455 }456 CATCH457 458 /*!459 Compute mapping of global index of grid source and grid destination460 \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1.461 \param [in] srcRank rank of client from which we demand global index of element source462 \param [in] src2DstMap mapping of global index of element source and global index of element destination463 \param [in] gridSrc Grid source464 \param [in] gridDst Grid destination465 \param [in] globalElementIndexOnProc Global index of element source on different client rank466 \param [out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination467 */468 void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid,469 const std::vector<int>& srcRank,470 std::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap,471 CGrid* gridSrc,472 CGrid* gridDst,473 std::vector<std::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc,474 SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)475 TRY476 {477 SourceDestinationIndexMap globaIndexWeightFromSrcToDst_tmp ;478 479 CContext* context = CContext::getCurrent();480 int clientRank = context->intraCommRank_;481 482 std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();483 std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();484 std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();485 CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;486 487 size_t nbElement = axisDomainSrcOrder.numElements();488 std::vector<size_t> nGlobSrc(nbElement);489 size_t globalSrcSize = 1;490 int domainIndex = 0, axisIndex = 0, scalarIndex = 0;491 for (int idx = 0; idx < nbElement; ++idx)492 {493 nGlobSrc[idx] = globalSrcSize;494 int elementDimension = axisDomainSrcOrder(idx);495 496 // If this is a domain497 if (2 == elementDimension)498 {499 globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();500 ++domainIndex;501 }502 else if (1 == elementDimension) // So it's an axis503 {504 globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();505 ++axisIndex;506 }507 else508 {509 globalSrcSize *= 1;510 ++scalarIndex;511 }512 }513 514 std::vector<CDomain*> domainListDestP = gridDst->getDomains();515 std::vector<CAxis*> axisListDestP = gridDst->getAxis();516 std::vector<CScalar*> scalarListDestP = gridDst->getScalars();517 CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;518 519 std::vector<size_t> nGlobDst(nbElement);520 size_t globalDstSize = 1;521 domainIndex = axisIndex = scalarIndex = 0;522 set<size_t> globalSrcStoreIndex ;523 524 for (int idx = 0; idx < nbElement; ++idx)525 {526 nGlobDst[idx] = globalDstSize;527 int elementDimension = axisDomainDstOrder(idx);528 529 // If this is a domain530 if (2 == elementDimension)531 {532 globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();533 ++domainIndex;534 }535 else if (1 == elementDimension) // So it's an axis536 {537 globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();538 ++axisIndex;539 }540 else541 {542 globalDstSize *= 1;543 ++scalarIndex;544 }545 }546 547 std::map< std::pair<size_t,size_t>, int > rankMap ;548 std::map< std::pair<size_t,size_t>, int >:: iterator rankMapIt ;549 550 for (int i = 0; i < srcRank.size(); ++i)551 {552 size_t ssize = 1;553 int rankSrc = srcRank[i];554 for (int idx = 0; idx < nbElement; ++idx)555 {556 ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();557 }558 559 std::vector<int> idxLoop(nbElement,0);560 std::vector<int> currentIndexSrc(nbElement, 0);561 std::vector<int> currentIndexDst(nbElement, 0);562 int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();563 size_t idx = 0;564 while (idx < ssize)565 {566 for (int ind = 0; ind < nbElement; ++ind)567 {568 if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())569 {570 idxLoop[ind] = 0;571 ++idxLoop[ind+1];572 }573 574 currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];575 }576 577 for (int ind = 0; ind < innnerLoopSize; ++ind)578 {579 currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];580 int globalElementDstIndexSize = 0;581 if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))582 {583 globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();584 }585 586 std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);587 size_t globalSrcIndex = 0;588 for (int idxElement = 0; idxElement < nbElement; ++idxElement)589 {590 if (idxElement == elementPositionInGrid)591 {592 for (int k = 0; k < globalElementDstIndexSize; ++k)593 {594 globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];595 }596 }597 else598 {599 for (int k = 0; k < globalElementDstIndexSize; ++k)600 {601 globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];602 }603 }604 globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];605 }606 607 for (int k = 0; k < globalElementDstIndexSize; ++k)608 {609 610 globaIndexWeightFromSrcToDst_tmp[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));611 rankMapIt=rankMap.find(make_pair(globalSrcIndex,globalDstVecIndex[k])) ;612 if (rankMapIt==rankMap.end()) rankMap[make_pair(globalSrcIndex,globalDstVecIndex[k])] = rankSrc ;613 else if (rankSrc==clientRank) rankMapIt->second = rankSrc ;614 }615 ++idxLoop[0];616 }617 idx += innnerLoopSize;618 }619 }620 621 // eliminate redondant global src point owned by differrent processes.622 // Avoid as possible to tranfer data from an other process if the src point is also owned by current process623 int rankSrc ;624 size_t globalSrcIndex ;625 size_t globalDstIndex ;626 double weight ;627 628 SourceDestinationIndexMap::iterator rankIt,rankIte ;629 std::unordered_map<size_t, std::vector<std::pair<size_t,double> > >::iterator globalSrcIndexIt, globalSrcIndexIte ;630 std::vector<std::pair<size_t,double> >::iterator vectIt,vectIte ;631 632 rankIt=globaIndexWeightFromSrcToDst_tmp.begin() ; rankIte=globaIndexWeightFromSrcToDst_tmp.end() ;633 for(;rankIt!=rankIte;++rankIt)634 {635 rankSrc = rankIt->first ;636 globalSrcIndexIt = rankIt->second.begin() ; globalSrcIndexIte = rankIt->second.end() ;637 for(;globalSrcIndexIt!=globalSrcIndexIte;++globalSrcIndexIt)638 {639 globalSrcIndex = globalSrcIndexIt->first ;640 vectIt = globalSrcIndexIt->second.begin() ; vectIte = globalSrcIndexIt->second.end() ;641 for(vectIt; vectIt!=vectIte; vectIt++)642 {643 globalDstIndex = vectIt->first ;644 weight = vectIt->second ;645 if (eliminateRedondantSrc_)646 {647 if (rankMap[make_pair(globalSrcIndex,globalDstIndex)] == rankSrc)648 globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ;649 }650 else globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ;651 }652 }653 }654 }655 CATCH656 657 /*!658 Find out proc and global index of axis source which axis destination is on demande659 \param[in] scalar Scalar destination660 \param[in] scalar Scalar source661 \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.662 \param[out] globalScalarIndexOnProc Global index of axis source on different procs663 */664 void CGenericAlgorithmTransformation::computeExchangeScalarIndex(CScalar* scalarDst,665 CScalar* scalarSrc,666 CArray<size_t,1>& destGlobalIndexPositionInGrid,667 std::unordered_map<int,std::vector<size_t> >& globalScalarIndexOnProc)668 TRY669 {670 CContext* context = CContext::getCurrent();671 int clientRank = context->intraCommRank_;672 int clientSize = context->intraCommSize_;673 674 globalScalarIndexOnProc.rehash(std::ceil(clientSize/globalScalarIndexOnProc.max_load_factor()));675 for (int idx = 0; idx < clientSize; ++idx)676 {677 globalScalarIndexOnProc[idx].push_back(0);678 }679 }680 CATCH681 682 /*!683 Find out proc and global index of axis source which axis destination is on demande684 \param[in] axisDst Axis destination685 \param[in] axisSrc Axis source686 \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.687 \param[out] globalAxisIndexOnProc Global index of axis source on different procs688 */689 void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,690 CAxis* axisSrc,691 CArray<size_t,1>& destGlobalIndexPositionInGrid,692 std::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)693 TRY694 {695 CContext* context = CContext::getCurrent();696 int clientRank = context->intraCommRank_;697 int clientSize = context->intraCommSize_;698 699 size_t globalIndex;700 int nIndexSize = axisSrc->index.numElements();701 CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;702 globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));703 for (int idx = 0; idx < nIndexSize; ++idx)704 {705 if (axisSrc->mask(idx))706 {707 globalIndex = axisSrc->index(idx);708 globalIndex2ProcRank[globalIndex].push_back(clientRank);709 }710 }711 712 CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);713 CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());714 for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)715 {716 globalAxisIndex(idx) = axisDst->index(idx);717 }718 dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);719 720 std::vector<int> countIndex(clientSize,0);721 const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();722 CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,723 ite = computedGlobalIndexOnProc.end();724 for (it = itb; it != ite; ++it)725 {726 const std::vector<int>& procList = it->second;727 for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];728 }729 730 globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));731 for (int idx = 0; idx < clientSize; ++idx)732 {733 if (0 != countIndex[idx])734 {735 globalAxisIndexOnProc[idx].resize(countIndex[idx]);736 countIndex[idx] = 0;737 }738 }739 740 for (it = itb; it != ite; ++it)741 {742 const std::vector<int>& procList = it->second;743 for (int idx = 0; idx < procList.size(); ++idx)744 {745 globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;746 ++countIndex[procList[idx]];747 }748 }749 }750 CATCH751 752 /*!753 Find out proc and global index of domain source which domain destination is on demande754 \param[in] domainDst Domain destination755 \param[in] domainSrc Domain source756 \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.757 \param[out] globalDomainIndexOnProc Global index of domain source on different procs758 */759 void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,760 CDomain* domainSrc,761 CArray<size_t,1>& destGlobalIndexPositionInGrid,762 std::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)763 TRY764 {765 CContext* context = CContext::getCurrent();766 int clientRank = context->intraCommRank_;767 int clientSize = context->intraCommSize_;768 769 int niGlobSrc = domainSrc->ni_glo.getValue();770 size_t globalIndex;771 int i_ind, j_ind;772 int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()773 : destGlobalIndexPositionInGrid.numElements();774 CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;775 globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));776 777 if (destGlobalIndexPositionInGrid.isEmpty())778 {779 for (int idx = 0; idx < nIndexSize; ++idx)780 {781 i_ind=domainSrc->i_index(idx) ;782 j_ind=domainSrc->j_index(idx) ;783 784 if (domainSrc->localMask(idx))785 {786 globalIndex = i_ind + j_ind * niGlobSrc;787 globalIndex2ProcRank[globalIndex].resize(1);788 globalIndex2ProcRank[globalIndex][0] = clientRank;789 }790 }791 }792 else793 {794 for (int idx = 0; idx < nIndexSize; ++idx)795 {796 // if (domainSrc->localMask(idx)) -> not necessairy, mask seem to be included in destGlobalIndexPositionInGrid(idx) (ym)797 globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);798 }799 }800 801 CArray<size_t,1> globalDomainIndex;802 if (destGlobalIndexPositionInGrid.isEmpty())803 {804 int niGlobDst = domainDst->ni_glo.getValue();805 globalDomainIndex.resize(domainDst->i_index.numElements());806 nIndexSize = domainDst->i_index.numElements();807 808 for (int idx = 0; idx < nIndexSize; ++idx)809 {810 i_ind=domainDst->i_index(idx) ;811 j_ind=domainDst->j_index(idx) ;812 813 globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;814 }815 }816 else817 {818 globalDomainIndex.reference(destGlobalIndexPositionInGrid);819 }820 821 CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);822 dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);823 824 std::vector<int> countIndex(clientSize,0);825 const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();826 CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,827 ite = computedGlobalIndexOnProc.end();828 for (it = itb; it != ite; ++it)829 {830 const std::vector<int>& procList = it->second;831 for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];832 }833 834 globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));835 for (int idx = 0; idx < clientSize; ++idx)836 {837 if (0 != countIndex[idx])838 {839 globalDomainIndexOnProc[idx].resize(countIndex[idx]);840 countIndex[idx] = 0;841 }842 }843 844 for (it = itb; it != ite; ++it)845 {846 const std::vector<int>& procList = it->second;847 for (int idx = 0; idx < procList.size(); ++idx)848 {849 globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;850 ++countIndex[procList[idx]];851 }852 }853 }854 CATCH855 856 void CGenericAlgorithmTransformation::computeTransformationMappingNonDistributed(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst,857 vector<int>& localSrc, vector<int>& localDst, vector<double>& weight,858 int& nlocalIndexDest)859 TRY860 {861 862 CContext* context = CContext::getCurrent();863 int nbClient = context->intraCommSize_;864 865 computePositionElements(gridDst, gridSrc);866 std::vector<CScalar*> scalarListDstP = gridDst->getScalars();867 std::vector<CAxis*> axisListDstP = gridDst->getAxis();868 std::vector<CDomain*> domainListDstP = gridDst->getDomains();869 CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;870 std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();871 std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();872 std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();873 CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;874 875 int nElement=axisDomainSrcOrder.numElements() ;876 int indSrc=1 ;877 int indDst=1 ;878 vector<int> nIndexSrc(nElement) ;879 vector<int> nIndexDst(nElement) ;880 vector< CArray<bool,1>* > maskSrc(nElement) ;881 vector< CArray<bool,1>* > maskDst(nElement) ;882 883 int nlocalIndexSrc=1 ;884 // int nlocalIndexDest=1 ;885 nlocalIndexDest=1 ;886 CArray<bool,1> maskScalar(1) ;887 maskScalar = true ;888 889 890 for(int i=0 ; i<nElement; i++)891 {892 int dimElement = axisDomainSrcOrder(i);893 if (2 == dimElement) //domain894 {895 CDomain* domain=domainListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ;896 nIndexSrc[i] = domain->i_index.numElements() ;897 maskSrc[i]=&domain->localMask ;898 }899 else if (1 == dimElement) //axis900 {901 CAxis* axis=axisListSrcP[elementPositionInGridSrc2AxisPosition_[i]] ;902 nIndexSrc[i] = axis->index.numElements() ;903 maskSrc[i]=&axis->mask ;904 }905 else //scalar906 {907 nIndexSrc[i]=1 ;908 maskSrc[i]=&maskScalar ;909 }910 nlocalIndexSrc=nlocalIndexSrc*nIndexSrc[i] ;911 }912 913 914 915 int offset=1 ;916 for(int i=0 ; i<nElement; i++)917 {918 int dimElement = axisDomainDstOrder(i);919 if (2 == dimElement) //domain920 {921 CDomain* domain=domainListDstP[elementPositionInGridDst2DomainPosition_[i]] ;922 int nIndex=domain->i_index.numElements() ;923 CArray<bool,1>& localMask=domain->localMask ;924 int nbInd=0 ;925 for(int j=0;j<nIndex;j++) if (localMask(j)) nbInd++ ;926 nIndexDst[i] = nbInd ;927 maskDst[i]=&domain->localMask ;928 }929 else if (1 == dimElement) //axis930 {931 CAxis* axis = axisListDstP[elementPositionInGridDst2AxisPosition_[i]] ;932 int nIndex=axis->index.numElements() ;933 CArray<bool,1>& localMask=axis->mask ;934 int nbInd=0 ;935 for(int j=0;j<nIndex;j++) if (localMask(j)) nbInd++ ;936 nIndexDst[i] = nbInd ;937 maskDst[i]=&axis->mask ;938 }939 else //scalar940 {941 nIndexDst[i]=1 ;942 maskDst[i]=&maskScalar ;943 }944 if (i<elementPositionInGrid) offset=offset*nIndexDst[i] ;945 nlocalIndexDest=nlocalIndexDest*nIndexDst[i] ;946 }947 948 vector<int> dstLocalInd ;949 int dimElement = axisDomainDstOrder(elementPositionInGrid);950 if (2 == dimElement) //domain951 {952 CDomain* domain = domainListDstP[elementPositionInGridDst2DomainPosition_[elementPositionInGrid]] ;953 int ni_glo=domain->ni_glo ;954 int nj_glo=domain->nj_glo ;955 int nindex_glo=ni_glo*nj_glo ;956 dstLocalInd.resize(nindex_glo,-1) ;957 int nIndex=domain->i_index.numElements() ;958 CArray<bool,1>& localMask=domain->localMask ;959 int unmaskedInd=0 ;960 int globIndex ;961 for(int i=0;i<nIndex;i++)962 {963 if (localMask(i))964 {965 globIndex=domain->j_index(i)*ni_glo+domain->i_index(i) ;966 dstLocalInd[globIndex]=unmaskedInd ;967 unmaskedInd++ ;968 }969 }970 }971 else if (1 == dimElement) //axis972 {973 CAxis* axis = axisListDstP[elementPositionInGridDst2AxisPosition_[elementPositionInGrid]] ;974 int nindex_glo=axis->n_glo ;975 dstLocalInd.resize(nindex_glo,-1) ;976 int nIndex=axis->index.numElements() ;977 CArray<bool,1>& localMask=axis->mask ; // axis mask must include later data_index978 int unmaskedInd=0 ;979 for(int i=0;i<nIndex;i++)980 {981 if (localMask(i))982 {983 dstLocalInd[axis->index(i)]=unmaskedInd ;984 unmaskedInd++ ;985 }986 }987 }988 else //scalar989 {990 dstLocalInd.resize(1) ;991 dstLocalInd[0]=0 ;992 }993 994 vector<vector<vector<pair<int,double> > > > dstIndWeight(transformationMapping_.size()) ;995 996 for(int t=0;t<transformationMapping_.size();++t)997 {998 TransformationIndexMap::const_iterator itTransMap = transformationMapping_[t].begin(),999 iteTransMap = transformationMapping_[t].end();1000 TransformationWeightMap::const_iterator itTransWeight = transformationWeight_[t].begin();1001 dstIndWeight[t].resize(nIndexSrc[elementPositionInGrid]) ;1002 1003 for(;itTransMap!=iteTransMap;++itTransMap,++itTransWeight)1004 {1005 int dst=dstLocalInd[itTransMap->first] ;1006 if (dst!=-1)1007 {1008 const vector<int>& srcs=itTransMap->second;1009 const vector<double>& weights=itTransWeight->second;1010 for(int i=0;i<srcs.size() ;i++) dstIndWeight[t][srcs[i]].push_back(pair<int,double>(dst*offset+t,weights[i])) ;1011 }1012 }1013 }1014 int srcInd=0 ;1015 int currentInd ;1016 int t=0 ;1017 int srcIndCompressed=0 ;1018 1019 nonDistributedrecursiveFunct(nElement-1,true,elementPositionInGrid,maskSrc,maskDst, srcInd, srcIndCompressed, nIndexSrc, t, dstIndWeight,1020 currentInd,localSrc,localDst,weight);1021 1022 }1023 CATCH1024 1025 1026 void CGenericAlgorithmTransformation::nonDistributedrecursiveFunct(int currentPos, bool masked, int elementPositionInGrid,1027 vector< CArray<bool,1>* >& maskSrc, vector< CArray<bool,1>* >& maskDst,1028 int& srcInd, int& srcIndCompressed, vector<int>& nIndexSrc,1029 int& t, vector<vector<vector<pair<int,double> > > >& dstIndWeight, int currentInd,1030 vector<int>& localSrc, vector<int>& localDst, vector<double>& weight)1031 TRY1032 {1033 int masked_ ;1034 if (currentPos!=elementPositionInGrid)1035 {1036 if (currentPos!=0)1037 {1038 CArray<bool,1>& mask = *maskSrc[currentPos] ;1039 1040 for(int i=0;i<nIndexSrc[currentPos];i++)1041 {1042 masked_=masked ;1043 if (!mask(i)) masked_=false ;1044 nonDistributedrecursiveFunct(currentPos-1, masked_, elementPositionInGrid, maskSrc, maskDst, srcInd, srcIndCompressed, nIndexSrc, t,1045 dstIndWeight, currentInd, localSrc, localDst, weight);1046 }1047 }1048 else1049 {1050 CArray<bool,1>& mask = *maskSrc[currentPos] ;1051 for(int i=0;i<nIndexSrc[currentPos];i++)1052 {1053 if (masked && mask(i))1054 {1055 if (dstIndWeight[t][currentInd].size()>0)1056 {1057 for(vector<pair<int,double> >::iterator it = dstIndWeight[t][currentInd].begin(); it!=dstIndWeight[t][currentInd].end(); ++it)1058 {1059 localSrc.push_back(srcIndCompressed) ;1060 localDst.push_back(it->first) ;1061 weight.push_back(it->second) ;1062 (it->first)++ ;1063 }1064 }1065 if (t < dstIndWeight.size()-1) t++ ;1066 srcIndCompressed ++ ;1067 }1068 srcInd++ ;1069 }1070 }1071 }1072 else1073 {1074 1075 if (currentPos!=0)1076 {1077 1078 CArray<bool,1>& mask = *maskSrc[currentPos] ;1079 for(int i=0;i<nIndexSrc[currentPos];i++)1080 {1081 t=0 ;1082 masked_=masked ;1083 if (!mask(i)) masked_=false ;1084 nonDistributedrecursiveFunct(currentPos-1, masked_, elementPositionInGrid, maskSrc, maskDst, srcInd,1085 srcIndCompressed, nIndexSrc, t, dstIndWeight , i, localSrc, localDst, weight);1086 }1087 }1088 else1089 {1090 for(int i=0;i<nIndexSrc[currentPos];i++)1091 {1092 if (masked)1093 {1094 t=0 ;1095 if (dstIndWeight[t][i].size()>0)1096 {1097 for(vector<pair<int,double> >::iterator it = dstIndWeight[t][i].begin(); it!=dstIndWeight[t][i].end(); ++it)1098 {1099 localSrc.push_back(srcIndCompressed) ;1100 localDst.push_back(it->first) ;1101 weight.push_back(it->second) ;1102 (it->first)++ ;1103 }1104 }1105 if (t < dstIndWeight.size()-1) t++ ;1106 srcIndCompressed ++ ;1107 }1108 srcInd++ ;1109 }1110 }1111 }1112 1113 }1114 CATCH1115 1116 /*!1117 Compute index mapping between element source and element destination with an auxiliary inputs which determine1118 position of each mapped index in global index of grid destination.1119 \param [in] dataAuxInputs auxiliary inputs1120 */1121 void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)1122 TRY1123 {1124 computeIndexSourceMapping_(dataAuxInputs);1125 }1126 CATCH1127 1128 std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()1129 TRY1130 {1131 return idAuxInputs_;1132 }1133 CATCH1134 1135 CGenericAlgorithmTransformation::AlgoTransType CGenericAlgorithmTransformation::type()1136 TRY1137 {1138 return type_;1139 }1140 CATCH1141 29 1142 30 … … 1154 42 CTransformFilter* CGenericAlgorithmTransformation::createTransformFilter(CGarbageCollector& gc, CGridAlgorithm* algo, bool detectMissingValues, double defaultValue) 1155 43 { 1156 return new CTransformFilter(gc, algo, detectMissingValues, defaultValue) ;44 return new CTransformFilter(gc, 1, algo, detectMissingValues, defaultValue) ; 1157 45 } 1158 46 1159 v oid CGenericAlgorithmTransformation::apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)47 vector<string> CGenericAlgorithmTransformation::getAuxFieldId(void) 1160 48 { 1161 // tranform into pure virtual funtion later 1162 abort() ; 49 return vector<string>() ; 1163 50 } 51 1164 52 }
Note: See TracChangeset
for help on using the changeset viewer.