#include "spatial_transform_filter.hpp" #include "grid_transformation.hpp" #include "context.hpp" #include "context_client.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, boost::shared_ptr > CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid, bool hasMissingValue, double missingValue) { if (!srcGrid || !destGrid) ERROR("std::pair, boost::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."); boost::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; boost::shared_ptr filter(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); } 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()) { boost::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) { CContextClient* client = CContext::getCurrent()->client; // 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::list >& listLocalIndexMaskOnDest = gridTransformation->getLocalMaskIndexOnGridDest(); 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::list >::const_iterator itLocalMaskIndexOnDest = listLocalIndexMaskOnDest.begin(); std::vector::const_iterator itAlgo = listAlgos.begin(); for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itLocalMaskIndexOnDest, ++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()); for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) { if (0 != itSend->second.numElements()) sendBuff[idxSendBuff] = 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(); 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()); } // 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) 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; 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); const std::vector& localMaskDest = *itLocalMaskIndexOnDest; for (int i = 0; i < localMaskDest.size(); ++i) if (localMaskDest[i]) dataCurrentDest(i) = 0.0; else dataCurrentDest(i) = defaultValue; std::vector localInitFlag(dataCurrentDest.numElements(), true); currentBuff = 0; bool firstPass=true; for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) { int countSize = itRecv->second.size(); const std::vector >& localIndex_p = itRecv->second; (*itAlgo)->apply(localIndex_p, recvBuff+currentBuff, dataCurrentDest, localInitFlag, ignoreMissingValue,firstPass); currentBuff += countSize; firstPass=false ; } (*itAlgo)->updateData(dataCurrentDest); idxSendBuff = 0; for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff) { if (0 != itSend->second.numElements()) delete [] sendBuff[idxSendBuff]; } 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; } } // namespace xios