#include "spatial_transform_filter.hpp" #include "grid_transformation.hpp" #include "context.hpp" #include "context_client.hpp" #include "timer.hpp" namespace xios { CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, double outputValue, size_t inputSlotsCount) : CFilter(gc, inputSlotsCount, engine), outputDefaultValue(outputValue) { /* Nothing to do */ } std::pair, std::shared_ptr > CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid, bool hasMissingValue, double missingValue) { if (!srcGrid || !destGrid) ERROR("std::pair, std::shared_ptr >" "buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)", "Impossible to build the filter graph if either the source or the destination grid are null."); std::shared_ptr firstFilter, lastFilter; // Note that this loop goes from the last transformation to the first transformation do { CGridTransformation* gridTransformation = destGrid->getTransformations(); CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations()); const std::vector& auxInputs = gridTransformation->getAuxInputs(); size_t inputCount = 1 + (auxInputs.empty() ? 0 : auxInputs.size()); double defaultValue = (hasMissingValue) ? std::numeric_limits::quiet_NaN() : 0.0; const CGridTransformationSelector::ListAlgoType& algoList = gridTransformation->getAlgoList() ; CGridTransformationSelector::ListAlgoType::const_iterator it ; bool isSpatialTemporal=false ; for (it=algoList.begin();it!=algoList.end();++it) if (it->second.first == TRANS_TEMPORAL_SPLITTING) isSpatialTemporal=true ; std::shared_ptr filter ; if( isSpatialTemporal) filter = std::shared_ptr(new CSpatialTemporalFilter(gc, engine, gridTransformation, defaultValue, inputCount)); else filter = std::shared_ptr(new CSpatialTransformFilter(gc, engine, defaultValue, inputCount)); if (!lastFilter) lastFilter = filter; else filter->connectOutput(firstFilter, 0); firstFilter = filter; for (size_t idx = 0; idx < auxInputs.size(); ++idx) { CField* fieldAuxInput = CField::get(auxInputs[idx]); fieldAuxInput->buildFilterGraph(gc, false); fieldAuxInput->getInstantDataFilter()->connectOutput(firstFilter,idx+1); } destGrid = gridTransformation->getGridSource(); } while (destGrid != srcGrid); return std::make_pair(firstFilter, lastFilter); } void CSpatialTransformFilter::onInputReady(std::vector data) { CSpatialTransformFilterEngine* spaceFilter = static_cast(engine); CDataPacketPtr outputPacket = spaceFilter->applyFilter(data, outputDefaultValue); if (outputPacket) onOutputReady(outputPacket); } CSpatialTemporalFilter::CSpatialTemporalFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, CGridTransformation* gridTransformation, double outputValue, size_t inputSlotsCount) : CSpatialTransformFilter(gc, engine, outputValue, inputSlotsCount), record(0) { const CGridTransformationSelector::ListAlgoType& algoList = gridTransformation->getAlgoList() ; CGridTransformationSelector::ListAlgoType::const_iterator it ; int pos ; for (it=algoList.begin();it!=algoList.end();++it) if (it->second.first == TRANS_TEMPORAL_SPLITTING) { pos=it->first ; if (pos < algoList.size()-1) ERROR("SpatialTemporalFilter::CSpatialTemporalFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, CGridTransformation* gridTransformation, double outputValue, size_t inputSlotsCount))", "temporal splitting operation must be the last of whole transformation on same grid") ; } CGrid* grid=gridTransformation->getGridDestination() ; CAxis* axis = grid->getAxis(gridTransformation->getElementPositionInGridDst2AxisPosition().find(pos)->second) ; nrecords = axis->index.numElements() ; } void CSpatialTemporalFilter::onInputReady(std::vector data) { CSpatialTransformFilterEngine* spaceFilter = static_cast(engine); CDataPacketPtr outputPacket = spaceFilter->applyFilter(data, outputDefaultValue); if (outputPacket) { size_t nelements=outputPacket->data.numElements() ; if (!tmpData.numElements()) { tmpData.resize(nelements); tmpData=outputDefaultValue ; } nelements/=nrecords ; size_t offset=nelements*record ; for(size_t i=0;idata(i) ; record ++ ; if (record==nrecords) { record=0 ; CDataPacketPtr packet = CDataPacketPtr(new CDataPacket); packet->date = data[0]->date; packet->timestamp = data[0]->timestamp; packet->status = data[0]->status; packet->data.resize(tmpData.numElements()); packet->data = tmpData; onOutputReady(packet); tmpData.resize(0) ; } } } CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation) : gridTransformation(gridTransformation) { if (!gridTransformation) ERROR("CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)", "Impossible to construct a spatial transform filter engine without a valid grid transformation."); } std::map > CSpatialTransformFilterEngine::engines; CSpatialTransformFilterEngine* CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation) { if (!gridTransformation) ERROR("CSpatialTransformFilterEngine& CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)", "Impossible to get the requested engine, the grid transformation is invalid."); std::map >::iterator it = engines.find(gridTransformation); if (it == engines.end()) { std::shared_ptr engine(new CSpatialTransformFilterEngine(gridTransformation)); it = engines.insert(std::make_pair(gridTransformation, engine)).first; } return it->second.get(); } CDataPacketPtr CSpatialTransformFilterEngine::apply(std::vector data) { /* Nothing to do */ } CDataPacketPtr CSpatialTransformFilterEngine::applyFilter(std::vector data, double defaultValue) { CDataPacketPtr packet(new CDataPacket); packet->date = data[0]->date; packet->timestamp = data[0]->timestamp; packet->status = data[0]->status; if (packet->status == CDataPacket::NO_ERROR) { if (1 < data.size()) // Dynamical transformations { std::vector* > dataAuxInputs(data.size()-1); for (size_t idx = 0; idx < dataAuxInputs.size(); ++idx) dataAuxInputs[idx] = &(data[idx+1]->data); gridTransformation->computeAll(dataAuxInputs, packet->timestamp); } packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements()); if (0 != packet->data.numElements()) (packet->data)(0) = defaultValue; apply(data[0]->data, packet->data); } return packet; } void CSpatialTransformFilterEngine::apply(const CArray& dataSrc, CArray& dataDest) { CTimer::get("CSpatialTransformFilterEngine::apply").resume(); CContextClient* client = CContext::getCurrent()->client; int rank; MPI_Comm_rank (client->intraComm, &rank); // Get default value for output data bool ignoreMissingValue = false; double defaultValue = std::numeric_limits::quiet_NaN(); if (0 != dataDest.numElements()) ignoreMissingValue = NumTraits::isNan(dataDest(0)); const std::list& listLocalIndexSend = gridTransformation->getLocalIndexToSendFromGridSource(); const std::list& listLocalIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest(); const std::list& listNbLocalIndexToReceive = gridTransformation->getNbLocalIndexToReceiveOnGridDest(); const std::vector& listAlgos = gridTransformation->getAlgos(); CArray dataCurrentDest(dataSrc.copy()); std::list::const_iterator itListSend = listLocalIndexSend.begin(), iteListSend = listLocalIndexSend.end(); std::list::const_iterator itListRecv = listLocalIndexToReceive.begin(); std::list::const_iterator itNbListRecv = listNbLocalIndexToReceive.begin(); std::vector::const_iterator itAlgo = listAlgos.begin(); for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itAlgo) { CArray dataCurrentSrc(dataCurrentDest); const CGridTransformation::SendingIndexGridSourceMap& localIndexToSend = *itListSend; // Sending data from field sources to do transformations std::map >::const_iterator itbSend = localIndexToSend.begin(), itSend, iteSend = localIndexToSend.end(); int idxSendBuff = 0; std::vector sendBuff(localIndexToSend.size()); double* sendBuffRank; for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) { int destRank = itSend->first; if (0 != itSend->second.numElements()) { if (rank != itSend->first) sendBuff[idxSendBuff] = new double[itSend->second.numElements()]; else sendBuffRank = new double[itSend->second.numElements()]; } } idxSendBuff = 0; std::vector sendRecvRequest; for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) { int destRank = itSend->first; const CArray& localIndex_p = itSend->second; int countSize = localIndex_p.numElements(); if (destRank != rank) { for (int idx = 0; idx < countSize; ++idx) { sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx)); } sendRecvRequest.push_back(MPI_Request()); MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest.back()); } else { for (int idx = 0; idx < countSize; ++idx) { sendBuffRank[idx] = dataCurrentSrc(localIndex_p(idx)); } } } // Receiving data on destination fields const CGridTransformation::RecvIndexGridDestinationMap& localIndexToReceive = *itListRecv; CGridTransformation::RecvIndexGridDestinationMap::const_iterator itbRecv = localIndexToReceive.begin(), itRecv, iteRecv = localIndexToReceive.end(); int recvBuffSize = 0; for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) { if (itRecv->first != rank ) recvBuffSize += itRecv->second.size(); } //(recvBuffSize < itRecv->second.size()) ? itRecv->second.size() : recvBuffSize; double* recvBuff; if (0 != recvBuffSize) recvBuff = new double[recvBuffSize]; int currentBuff = 0; for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) { int srcRank = itRecv->first; if (srcRank != rank) { int countSize = itRecv->second.size(); sendRecvRequest.push_back(MPI_Request()); MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest.back()); currentBuff += countSize; } } std::vector status(sendRecvRequest.size()); MPI_Waitall(sendRecvRequest.size(), &sendRecvRequest[0], &status[0]); dataCurrentDest.resize(*itNbListRecv); dataCurrentDest = 0.0; std::vector localInitFlag(dataCurrentDest.numElements(), true); currentBuff = 0; bool firstPass=true; for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) { const std::vector >& localIndex_p = itRecv->second; int srcRank = itRecv->first; if (srcRank != rank) { int countSize = itRecv->second.size(); (*itAlgo)->apply(localIndex_p, recvBuff+currentBuff, dataCurrentDest, localInitFlag, ignoreMissingValue,firstPass); currentBuff += countSize; } else { (*itAlgo)->apply(localIndex_p, sendBuffRank, dataCurrentDest, localInitFlag, ignoreMissingValue,firstPass); } firstPass=false ; } (*itAlgo)->updateData(dataCurrentDest); idxSendBuff = 0; for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) { if (0 != itSend->second.numElements()) { if (rank != itSend->first) delete [] sendBuff[idxSendBuff]; else delete [] sendBuffRank; } } if (0 != recvBuffSize) delete [] recvBuff; } if (dataCurrentDest.numElements() != dataDest.numElements()) ERROR("CSpatialTransformFilterEngine::apply(const CArray& dataSrc, CArray& dataDest)", "Incoherent between the received size and expected size. " << std::endl << "Expected size: " << dataDest.numElements() << std::endl << "Received size: " << dataCurrentDest.numElements()); dataDest = dataCurrentDest; CTimer::get("CSpatialTransformFilterEngine::apply").suspend() ; } } // namespace xios