source: XIOS/trunk/src/filter/spatial_transform_filter.cpp @ 842

Last change on this file since 842 was 842, checked in by mhnguyen, 5 years ago

Small improvements on transformation mapping

+) Remove complex structure of transformation mapping by a simpler one

Test
+) On Curie
+) All tests pass

File size: 9.0 KB
Line 
1#include "spatial_transform_filter.hpp"
2#include "grid_transformation.hpp"
3#include "context.hpp"
4#include "context_client.hpp"
5
6namespace xios
7{
8  CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, size_t inputSlotsCount)
9    : CFilter(gc, inputSlotsCount, engine)
10  { /* Nothing to do */ }
11
12  std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >
13  CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)
14  {
15    if (!srcGrid || !destGrid)
16      ERROR("std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >"
17            "buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)",
18            "Impossible to build the filter graph if either the source or the destination grid are null.");
19
20    boost::shared_ptr<CSpatialTransformFilter> firstFilter, lastFilter;
21    // Note that this loop goes from the last transformation to the first transformation
22    do
23    {
24      CGridTransformation* gridTransformation = destGrid->getTransformations();
25      CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations());
26      const std::vector<StdString>& auxInputs = gridTransformation->getAuxInputs();
27      size_t inputCount = 1 + (auxInputs.empty() ? 0 : auxInputs.size());
28      boost::shared_ptr<CSpatialTransformFilter> filter(new CSpatialTransformFilter(gc, engine, inputCount));
29
30      if (!lastFilter)
31        lastFilter = filter;
32      else
33        filter->connectOutput(firstFilter, 0);
34
35      firstFilter = filter;
36      for (size_t idx = 0; idx < auxInputs.size(); ++idx)
37      {
38        CField* fieldAuxInput = CField::get(auxInputs[idx]);
39        fieldAuxInput->buildFilterGraph(gc, false);
40        fieldAuxInput->getInstantDataFilter()->connectOutput(firstFilter,idx+1);
41      }
42
43      destGrid = gridTransformation->getGridSource();
44    }
45    while (destGrid != srcGrid);
46
47    return std::make_pair(firstFilter, lastFilter);
48  }
49
50  CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)
51    : gridTransformation(gridTransformation)
52  {
53    if (!gridTransformation)
54      ERROR("CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)",
55            "Impossible to construct a spatial transform filter engine without a valid grid transformation.");
56  }
57
58  std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > CSpatialTransformFilterEngine::engines;
59
60  CSpatialTransformFilterEngine* CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)
61  {
62    if (!gridTransformation)
63      ERROR("CSpatialTransformFilterEngine& CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)",
64            "Impossible to get the requested engine, the grid transformation is invalid.");
65
66    std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines.find(gridTransformation);
67    if (it == engines.end())
68    {
69      boost::shared_ptr<CSpatialTransformFilterEngine> engine(new CSpatialTransformFilterEngine(gridTransformation));
70      it = engines.insert(std::make_pair(gridTransformation, engine)).first;
71    }
72
73    return it->second.get();
74  }
75
76  CDataPacketPtr CSpatialTransformFilterEngine::apply(std::vector<CDataPacketPtr> data)
77  {
78    CDataPacketPtr packet(new CDataPacket);
79    packet->date = data[0]->date;
80    packet->timestamp = data[0]->timestamp;
81    packet->status = data[0]->status;
82
83    if (packet->status == CDataPacket::NO_ERROR)
84    {
85      if (1 < data.size())  // Dynamical transformations
86      {
87        std::vector<CArray<double,1>* > dataAuxInputs(data.size()-1);
88        for (size_t idx = 0; idx < dataAuxInputs.size(); ++idx) dataAuxInputs[idx] = &(data[idx+1]->data);
89        gridTransformation->computeAll(dataAuxInputs, packet->timestamp);
90      }
91      packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements());
92      apply(data[0]->data, packet->data);
93    }
94
95    return packet;
96  }
97
98  void CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)
99  {
100    CContextClient* client = CContext::getCurrent()->client;
101
102    const std::list<CGridTransformation::SendingIndexGridSourceMap>& listLocalIndexSend = gridTransformation->getLocalIndexToSendFromGridSource();
103    const std::list<CGridTransformation::RecvIndexGridDestinationMap>& listLocalIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest();
104    const std::list<size_t>& listNbLocalIndexToReceive = gridTransformation->getNbLocalIndexToReceiveOnGridDest();
105
106    CArray<double,1> dataCurrentDest(dataSrc.copy());
107
108    std::list<CGridTransformation::SendingIndexGridSourceMap>::const_iterator itListSend  = listLocalIndexSend.begin(),
109                                                                              iteListSend = listLocalIndexSend.end();
110    std::list<CGridTransformation::RecvIndexGridDestinationMap>::const_iterator itListRecv = listLocalIndexToReceive.begin();
111    std::list<size_t>::const_iterator itNbListRecv = listNbLocalIndexToReceive.begin();
112
113    for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv)
114    {
115      CArray<double,1> dataCurrentSrc(dataCurrentDest);
116      const CGridTransformation::SendingIndexGridSourceMap& localIndexToSend = *itListSend;
117
118      // Sending data from field sources to do transformations
119      std::map<int, CArray<int,1> >::const_iterator itbSend = localIndexToSend.begin(), itSend,
120                                                    iteSend = localIndexToSend.end();
121      int idxSendBuff = 0;
122      std::vector<double*> sendBuff(localIndexToSend.size());
123      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
124      {
125        if (0 != itSend->second.numElements())
126          sendBuff[idxSendBuff] = new double[itSend->second.numElements()];
127      }
128
129      idxSendBuff = 0;
130      std::vector<MPI_Request> sendRecvRequest;
131      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
132      {
133        int destRank = itSend->first;
134        const CArray<int,1>& localIndex_p = itSend->second;
135        int countSize = localIndex_p.numElements();
136        for (int idx = 0; idx < countSize; ++idx)
137        {
138          sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx));
139        }
140        sendRecvRequest.push_back(MPI_Request());
141        MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest.back());
142      }
143
144      // Receiving data on destination fields
145      const CGridTransformation::RecvIndexGridDestinationMap& localIndexToReceive = *itListRecv;
146      CGridTransformation::RecvIndexGridDestinationMap::const_iterator itbRecv = localIndexToReceive.begin(), itRecv,
147                                                                       iteRecv = localIndexToReceive.end();
148      int recvBuffSize = 0;
149      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += itRecv->second.size(); //(recvBuffSize < itRecv->second.size())
150                                                                       //? itRecv->second.size() : recvBuffSize;
151      double* recvBuff;
152      if (0 != recvBuffSize) recvBuff = new double[recvBuffSize];
153      int currentBuff = 0;
154      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
155      {
156        int srcRank = itRecv->first;
157        int countSize = itRecv->second.size();
158        sendRecvRequest.push_back(MPI_Request());
159        MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest.back());
160        currentBuff += countSize;
161      }
162      std::vector<MPI_Status> status(sendRecvRequest.size());
163      MPI_Waitall(sendRecvRequest.size(), &sendRecvRequest[0], &status[0]);
164
165      dataCurrentDest.resize(*itNbListRecv);
166      dataCurrentDest = 0.0;
167      currentBuff = 0;
168      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
169      {
170        int countSize = itRecv->second.size();
171        const std::vector<std::pair<int,double> >& localIndex_p = itRecv->second;
172        for (int idx = 0; idx < countSize; ++idx)
173        {
174          dataCurrentDest(localIndex_p[idx].first) += *(recvBuff+currentBuff+idx) * localIndex_p[idx].second;
175        }
176        currentBuff += countSize;
177      }
178
179      idxSendBuff = 0;
180      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
181      {
182        if (0 != itSend->second.numElements())
183          delete [] sendBuff[idxSendBuff];
184      }
185      if (0 != recvBuffSize) delete [] recvBuff;
186    }
187    if (dataCurrentDest.numElements() != dataDest.numElements())
188    ERROR("CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)",
189          "Incoherent between the received size and expected size" <<
190          "Expected size: " << dataDest.numElements() <<
191          "Received size: " << dataCurrentDest.numElements());
192
193    dataDest = dataCurrentDest;
194  }
195} // namespace xios
Note: See TracBrowser for help on using the repository browser.